evaluate.f90

! Part of Zinc FE package. Author: John Blackburn

module evaluate

! The user can assign values to parameters that can be used in expressions with 
! the subroutine defparam. The calling syntax is
!
!     call defparam(symbol,value) or call defparam(symbol,expr)
!
! where symbol is the desired parameter name; value is a real, integer, or 
! complex variable (single or double precision); and expr is a string
! containing an expression to be evaluated. The value obtained by evaluating the
! expression expr is associated with the parameter symbol. Parameter names must
! begin with a letter (a-z, A-Z) and must not be longer than 24 characters.
! Parameter names are not case dependent.
!
! An expression can be evaluated with the subroutine evalexpr. The calling
! syntax is 
!
!          call evalexpr(expr,value)
!
! where expr is a string containing the expression to be evaluated; value is the
! result (single or double precision real, complex or integer). The 
! expression can contain the arithmetic operations +, -, *, /, or ^ as well as
! the functions sin, cos, tan, log, ln, abs, exp, sqrt, real, imag, conjg, and 
! ang (the function ang calculates the phase angle of its complex argument). The
! expression can also contain numerical values and previously defined parameters
! Grouping by nested levels of parentheses is also allowed. The parameters pi
! and i (imaginary unit) are predefined. Complex numbers can be entered as a+i*b 
! if the parameter i has not been redefined by the user. Complex numbers can also
! be entered using complex(a,b).
! Example expression: 
!
!          conjg(((cos(x) + sqrt(a+i*b))^2+complex(ln(1.6e-4),20))/2)
!
! An equation of the form <symbol> = <expression> can be evaluated using the 
! subroutine evaleqn. The calling syntax is
!
!          call evaleqn(eqn)       
!
! where eqn is a string containing the equation. The right-hand-side of the
! equation is evaluated and assigned to the symbol given by the left-hand-side.
!
! The value assigned to a symbol can be retrieved using the subroutine getparam.
! The calling syntax is
!
!          call getparam(sym,value)
!
! where sym is a symbol string; value is a numeric variable (any of the six 
! standard types).
!
! The symbols and their values in the symbol table can be listed using the
! subroutine listvar. The variable ierr is always available following a call
! to any of the above subroutines and is zero if there were no errors. The
! possible nonzero values for ierr are
!
! 1       Expression empty
! 2       Parentheses don't match
! 3       Number string does not correspond to a valid number
! 4       Undefined symbol
! 5       Less than two operands for binary operation
! 6       No operand for unary plus or minus operators
! 7       No argument(s) for function
! 8       Zero or negative real argument for logarithm
! 9       Negative real argument for square root
! 10      Division by zero
! 11      Improper symbol format
! 12      Missing operator
! 13      Undefined function
! 14      Argument of tangent function a multiple of pi/2
! 15      Max no. of parameters exceeded [jfb]
! 16      Max no. of tokens exceeded [jfb]

use precision
use strings

implicit none

save
private
public :: valuep,evalexpr,defparam,evaleqn,getparam,listvar,evaltoken,tokenize,get_next_token
public :: ierr,nparams,params,numtok
public :: item,Tlist

type item                         
   character(len=24):: char
   character :: type
   complex(kr8) :: value
   integer :: ind
end type item

type param                        
  character (len=24):: symbol
  complex(kc8):: value
end type param

interface defparam                
  module procedure strdef       ! value given by expression      
  module procedure valdef_dc    ! Double precision complex value
  module procedure valdef_sc    ! Single precision complex value
  module procedure valdef_dr    ! Double precision real value
  module procedure valdef_sr    ! Single precision real value
  module procedure valdef_di    ! Double precision integer value
  module procedure valdef_si    ! Single precision integer value
end interface

interface evalexpr
  module procedure evalexpr_dc  ! Double precision complex result
  module procedure evalexpr_sc  ! Single precision complex result
  module procedure evalexpr_dr  ! Double precision real result
  module procedure evalexpr_sr  ! Single precision real result
  module procedure evalexpr_di  ! Double precision integer result
  module procedure evalexpr_si  ! Single precision integer result
end interface

interface getparam
  module procedure getparam_dc  ! Double precision complex result
  module procedure getparam_sc  ! Single precision complex result
  module procedure getparam_dr  ! Double precision real result
  module procedure getparam_sr  ! Single precision real result
  module procedure getparam_di  ! Double precision integer result
  module procedure getparam_si  ! Single precision integer result
end interface

interface evaltoken
   module procedure evaltoken_DC
   module procedure evaltoken_DR
end interface

integer, parameter :: numtok=350  ! Maximum number of tokens 
integer, parameter :: mxparams=1000 ! Maximum number of variables
type(param) :: params(mxparams)       ! Symbol table
integer :: nparams=0,itop,ibin
complex(kc8) :: valstack(numtok) ! Stack used in evaluation of expression
type(item):: opstack(0:numtok)     ! Operator stack used in conversion to postfix
integer :: ierr                  ! Error flag

type Tlist
   integer ntok
   type(item), allocatable :: token(:)
end type Tlist

!**********************************************************************

contains

!**********************************************************************


SUBROUTINE EVALEXPR_DC(expr,val)    ! Evaluate expression expr for
                                    ! val double precision complex

character (len=*),intent(in) :: expr
complex(kc8) :: val
character (len=len(expr)+1) :: tempstr
integer :: isp(0:numtok)          ! On stack priority of operators in opstack
integer :: lstr
complex(kc8) :: cval,oper1,oper2
real(kr8) :: valr,vali
type(item):: token(numtok)      ! List of tokens ( a token is an operator or 
                                ! operand) in postfix order
type(item) :: x,junk,tok
integer i,icp,insp,isum,ntok

ierr=0
token(1:)%char=' '

if(nparams == 0) then                  ! Initialize symbol table
  params(1)%symbol='PI'
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
  params(2)%symbol='I'
  params(2)%value=(0.0_kr8,1.0_kr8)
  nparams=2
end if

! ----------------------------------------------------------------------
! Check the expression does not have too many tokens [jfb]
! ----------------------------------------------------------------------

     tempstr=adjustl(expr)
     call removesp(tempstr)
     ntok=0
     do
        call get_next_token(tempstr,tok,icp,insp)
        if (tok%type.eq.'E') exit

        ntok=ntok+1
     enddo

     if (ntok.gt.numtok) then
        ierr=16
        write (*,*) 'Error: Too many tokens in expression'
        write (*,*) 'Found ',ntok,': maximum',numtok
        return
     endif

if(len_trim(expr) == 0) then           ! Expression empty
  ierr=1
  write(*,*) 'Error: expression being evaluated is empty'
  return
end if

tempstr=adjustl(expr)
call removesp(tempstr)   ! Removes spaces, tabs, and control characters

! ****************************************************************************
! STEP 1:  Convert string to token array. Each token is either an operator or
!          an operand. Token array will be in postfix (reverse Polish) order.
!*****************************************************************************

ntok=0
ibin=0
itop=0
isp(:)=0                 ! added by JFB
do
  lstr=len_trim(tempstr)
  call get_next_token(tempstr(1:lstr),tok,icp,insp)
  select case(tok%type)
  case('S')
    ntok=ntok+1
    token(ntok)=tok
!    print *,ntok,tok%char
  case('E')
    do 
      if(itop < 1)exit
      call popop(x)        ! Output remaining operators on stack
      ntok=ntok+1
      token(ntok)=x
    end do
    ntok=ntok+1
    token(ntok)=tok
    exit
  case('R')  ! Token is right parenenthesis
    do 
      if(opstack(itop)%type == 'L') exit  ! Output operators on stack down
      call popop(x)                       ! to left parenthesis
      ntok=ntok+1
      token(ntok)=x
    end do                             
    call popop(junk)                      ! Remove left parenthesis from stack
    if(opstack(itop)%type == 'F') then    ! Output function name if present
      call popop(x)
      ntok=ntok+1
      token(ntok)=x
    end if                               
  case('D')  ! Token is comma
    do 
      if(opstack(itop)%type == 'L') exit  ! Output operators on stack down
      call popop(x)                       ! to left parenthesis
      ntok=ntok+1
      token(ntok)=x
    end do                              
  case('U','B','L','F') ! Token is operator, left parenthesis or function name
    do 
      if(isp(itop) < icp) exit            ! Output operators on stack having
      call popop(x)                       ! an instack priority that is
      ntok=ntok+1                         ! greater than or equal to the
      token(ntok)=x                       ! priority of the incoming operator  
    end do                                
    call pushop(tok)     ! Put incoming operator on stack                       
    isp(itop)=insp
  end select
end do

isum=0                                 ! Error check for matching parentheses
do i=1,ntok
  if(token(i)%type == 'L' ) isum=isum+1
  if(token(i)%type == 'R' ) isum=isum-1
end do
if(isum /= 0) then
  ierr=2
  write(*,*) 'Error in the evaluation of the expression ',trim(expr)
  write(*,*) "Parentheses don't match"
  write(*,*)
  return
end if


!*****************************************************************************
! STEP 2: Evaluate token string in postfix order
!*****************************************************************************

itop=0
do i=1,ntok
  x=token(i)
  select case(x%type)
  case('E')  ! Token is end token
    if(itop>1) then                
      ierr=12
      write(*,*) 'Error: missing operator in expression ',trim(expr)
      write(*,*)
      return
    end if
    call popval(val)               ! Final result left on stack of values
    exit
  case('S')  ! Token is operand
    call valuep(x%char,cval)       ! Evaluate operand
    if(ierr/=0) return
    call pushval(cval)             ! Put value of operand on stack
  case('B')  ! Token is a binary operator
    if(itop < 2) then
      ierr=5
      write(*,*) 'Error in evaluation of expression ',trim(expr)
      write(*,*) 'Less than two operands for binary operator  '&
                 ,trim(x%char)
      write(*,*)
      return
    end if                         
    call popval(oper1)             ! Pull off top two values from stack
    call popval(oper2)
    select case(trim(x%char))      ! Perform operation on values
    case('^')
      cval=oper2**oper1
    case('*')
      cval=oper2*oper1
    case('/')
       if(oper1 == (0._kr8,0._kr8)) then
          ierr=10
          write(*,*) 'Warning in expression ',trim(expr)
          write(*,*) 'Division by zero'
          write(*,*)
          cval=0
       else
          cval=oper2/oper1
       end if
    case('+')
      cval=oper2+oper1
    case('-')
      cval=oper2-oper1
    end select
    call pushval(cval)             ! Put result back on stack
  case('U')  ! Token is unary operator
    if(itop == 0) then
      ierr=6
      write(*,*) 'Error in expression ',trim(expr)
      write(*,*) 'No operand for unary operator ',trim(x%char)
      write(*,*)
      return
    else
      call popval(oper1)           ! Pull top value off stack
    end if
    select case(trim(x%char))      ! Operate on value
    case('+')
      cval=oper1
    case('-')
      cval=-oper1
    end select
    call pushval(cval)             ! Put result back on stack
  case('F')  ! Token is a function name
    if(itop == 0) then
      ierr=7
      write(*,*) 'Error in expression ',trim(expr)
      write(*,*) 'Missing argument(s) for function ',trim(x%char)
      write(*,*)
      return
    else  
      call popval(oper1)           ! Pull top value off stack
    end if 
    tempstr=uppercase(x%char)
    select case(trim(tempstr))      ! Evaluate function
    case('SIN')
      cval=sin(oper1)
    case('COS')
      cval=cos(oper1)
    case('TAN')
      oper2=cos(oper1)
      if(abs(oper2) == 0.0_kr8) then
        ierr=14
        write(*,*) 'Error: argument of tan function a multiple',&
        ' of pi/2 in expression ',trim(expr)
        write(*,*)
        cval=0
      else 
        cval=sin(oper1)/oper2
      endif
    case('SQRT')
      if(real(oper1,kr8) < 0. .and. aimag(oper1)==0.) then
        ierr=9
        write(*,*) 'Warning: square root of negative real number',&
                   ' in expression ',trim(expr)
        write(*,*)
      end if
      cval=sqrt(oper1)
    case('ABS')
      cval=abs(oper1)
    case('LN')
      if(real(oper1,kr8) <= 0. .and. aimag(oper1)==0.) then
        ierr=8
        write(*,*) 'Error: negative real or zero argument for',&
                   ' natural logarithm in expression ',trim(expr)
        write(*,*)
        cval=0
      else
        cval=log(oper1)
      end if

    case('LOG')
      if(real(oper1,kr8) <= 0. .and. aimag(oper1)==0.) then
        ierr=8
        write(*,*) 'Error: negative real or zero argument for base',&
                   '10 logarithm in expression ',trim(expr)
        write(*,*)
        cval=0
      else
        cval=log(oper1)/2.30258509299405_kr8
      end if

    case('EXP')
      cval=exp(oper1)
    case('COMPLEX')
      if(itop == 0) then
        ierr=7
        write(*,*) 'Error in expression ',trim(expr)
        write(*,*) 'Missing argument(s) for function ',trim(x%char)
        write(*,*)
        return
      else  
        call popval(oper2)  ! Pull second argument off stack
      end if 
      valr=real(oper2,kr8)
      vali=real(oper1,kr8)
      cval=cmplx(valr,vali,kc8)
    case('CONJG')
      cval=conjg(oper1)
    case('ANG')
      cval=atan2(aimag(oper1),real(oper1,kr8))
    case('REAL')
      cval=real(oper1,kr8)
    case('IMAG')
      cval=aimag(oper1)
    case default ! Undefined function
      ierr=13
      write(*,*) 'Error: the function ',trim(x%char), ' is undefined',&
                 ' in the expression ',trim(expr)
      write(*,*)
      return
    end select
    call pushval(cval)    ! Put result back on stack
  end select
end do

end subroutine evalexpr_dc

!**********************************************************************

SUBROUTINE GET_NEXT_TOKEN(str,tok,icp,isp)

character(len=*) :: str
character :: cop,chtemp
type(item) :: tok
integer :: icp, isp
integer inext,ipos,lstr,ntok

lstr=len_trim(str)
if(lstr == 0) then
  tok%char='#'             ! Output end token
  tok%type='E'
  return
end if
ipos=scan(str,'+-*/^(),')  ! Look for an arithmetic operator 
                           ! + - * / ^ ( ) or ,
if (ipos.ne.0) cop=str(ipos:ipos)                 
select case (ipos)              
case(0)    ! Operators not present
  ntok=ntok+1
  tok%char=str
  tok%type='S'
  str=''
  icp=0
  isp=0
case(1) 
  tok%char=cop
  select case(cop)
  case('+','-')
    if(ibin==0) then
      tok%type='U'
      icp=4
      isp=3
    else
      tok%type='B'
      icp=1
      isp=1
    end if
    ibin=0
  case('*','/')
    tok%type='B'
    icp=2
    isp=2
    ibin=0
  case('^')
    tok%type='B'
    icp=4
    isp=3
    ibin=0
  case('(')
    tok%type='L'
    icp=4
    isp=0
    ibin=0
  case(')')
    tok%type='R'
    icp=0
    isp=0
    ibin=1
  case(',')
    tok%type='D'
    icp=0
    isp=0
    ibin=0
  end select
  str=str(2:)
case(2:)
  select case(cop)
  case('(')
    tok%char=str(1:ipos-1)
    tok%type='F'
    icp=4
    isp=0
    ibin=0
    str=str(ipos:)
  case('+','-')
    chtemp=uppercase(str(ipos-1:ipos-1))
    if(is_letter(str(1:1)) .or. chtemp/='E') then
      tok%char=str(1:ipos-1)
      tok%type='S'
      icp=0
      isp=0
      ibin=1
      str=str(ipos:)
    else
      inext=scan(str(ipos+1:),'+-*/^(),')
      if(inext==0) then
        tok%char=str
        tok%type='S'
        icp=0
        isp=0
        ibin=0
        str=''
      else
        tok%char=str(1:ipos+inext-1)
        tok%type='S'
        icp=0
        isp=0
        ibin=1
        str=str(ipos+inext:)
      end if
    end if
  case default
    tok%char=str(1:ipos-1)
    tok%type='S'
    icp=0
    isp=0
    ibin=1
    str=str(ipos:)
  end select
end select

end subroutine get_next_token


!**********************************************************************

SUBROUTINE EVALEXPR_SC(expr,val)    ! Evaluate expression expr for
                                    ! val single precision complex
character(len=*) :: expr
complex(kc4) :: val
complex(kc8) :: vald

call evalexpr_dc(expr,vald)
val=vald

end subroutine evalexpr_sc

!**********************************************************************

SUBROUTINE EVALEXPR_SR(expr,val)    ! Evaluate expression expr for
                                    ! val single precision real
character(len=*) :: expr
real(kr4) :: val
complex(kc8) :: vald

call evalexpr_dc(expr,vald)
val=real(vald)

end subroutine evalexpr_sr

!**********************************************************************

SUBROUTINE EVALEXPR_DR(expr,val)    ! Evaluate expression expr for 
                                    ! val double precision real
character(len=*) :: expr
real(kr8) :: val
complex(kc8) :: vald

call evalexpr_dc(expr,vald)
val=real(vald,kr8)

end subroutine evalexpr_dr

!**********************************************************************

SUBROUTINE EVALEXPR_SI(expr,ival)   ! Evaluate expression expr for 
                                    ! ival single precision integer
character(len=*) :: expr
integer(ki4) :: ival
complex(kc8) :: vald

call evalexpr_dc(expr,vald)
ival=nint(real(vald,kr8),ki4)

end subroutine evalexpr_si

!**********************************************************************

SUBROUTINE EVALEXPR_DI(expr,ival)   ! Evaluate expression expr for 
                                    ! ival double precision integer
character(len=*) :: expr
integer(ki8) :: ival
complex(kc8) :: vald

call evalexpr_dc(expr,vald)
ival=nint(real(vald,kr8),ki8)

end subroutine evalexpr_di


!**********************************************************************
SUBROUTINE VALDEF_DC(sym,val)    ! Associates sym with val in symbol table,
                                 ! val double precision complex
character(len=*) :: sym
character(len=len_trim(sym)) :: usym
complex(kc8) :: val

integer i

ierr=0
if(nparams == 0) then               ! Initialize symbol table
  params(1)%symbol='PI'
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
  params(2)%symbol='I'
  params(2)%value=(0.0_kr8,1.0_kr8)
  nparams=2
end if

! Assign val to sym if sym is already in symbol table
usym=uppercase(sym)
if((.not.is_letter(sym(1:1))) .or. len_trim(sym)>24) then
  ierr=11
  write(*,*) 'Error: symbol ',trim(sym),' has improper format'
  write(*,*)
  return
end if
do i=1,nparams
  if(trim(usym)==trim(params(i)%symbol)) then
    params(i)%value=val
    return
  end if
end do

if (nparams.ge.mxparams) then
   ierr=15
   write (*,*) 'Error: Too many variables'
   return
endif

nparams=nparams+1    ! Otherwise assign val to new symbol sym
params(nparams)%symbol=usym
params(nparams)%value=val

end subroutine valdef_dc


!**********************************************************************

SUBROUTINE VALDEF_SC(sym,val)     ! Associates sym with val in symbol table,
                                  ! val single precision complex
character(len=*) :: sym
complex(kc4) :: val
complex(kc8) :: vald

vald=val
call valdef_dc(sym,vald)

end subroutine valdef_sc


!**********************************************************************

SUBROUTINE VALDEF_DR(sym,val)    ! Associates sym with val in symbol table,
                                 ! val double precision real
character(len=*) :: sym
real(kr8) :: val
complex(kc8) :: vald

vald=cmplx(val,0.0_kr8,kc8)
call valdef_dc(sym,vald)

end subroutine valdef_dr


!**********************************************************************

SUBROUTINE VALDEF_SR(sym,val)    ! Associates sym with val in symbol table,
                                 ! val single precision real
character(len=*) :: sym
real(kr4) :: val
complex(kc8) :: vald

vald=cmplx(val,0.0,kc8)
call valdef_dc(sym,vald)

end subroutine valdef_sr


!**********************************************************************

SUBROUTINE VALDEF_DI(sym,ival)   ! Associates sym with ival in symbol table,
                                 ! ival double precision integer 
character(len=*) :: sym
integer(ki8) :: ival
complex(kc8) :: vald

vald=cmplx(real(ival,kr8),0.0_kr8,kc8)
call valdef_dc(sym,vald)

end subroutine valdef_di


!**********************************************************************

SUBROUTINE VALDEF_SI(sym,ival)   ! Associates sym with ival in symbol table,
                                 ! ival single precision integer
character(len=*) :: sym
integer(ki4) :: ival
complex(kc8) :: vald

vald=cmplx(real(ival,kr8),0.0,kc8)
call valdef_dc(sym,vald)

end subroutine valdef_si


!**********************************************************************

SUBROUTINE STRDEF(sym,expr)      ! Associates sym with the value of the
                                 ! expression expr

character(len=*) :: sym,expr
complex(kc8) :: val

if(nparams == 0) then            ! Initialize symbol table
  params(1)%symbol='PI'
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
  params(2)%symbol='I'
  params(2)%value=(0.0_kr8,1.0_kr8)
  nparams=2
end if

call evalexpr_dc(expr,val)       ! val is value of expression expr
if(ierr==0 .or. ierr==9) then
  call valdef_dc(sym,val)          ! Assign val to symbol sym
end if

end subroutine strdef


!**********************************************************************

SUBROUTINE VALUEP(xinchar,cval)  ! Finds double precision complex value 
                                 ! corresponding to number string xinchar 
                                 ! or value in symbol table corresponding 
                                 ! to symbol name xinchar.

character (len=*):: xinchar
complex(kc8) :: cval
real(kr8) :: rval
integer ios

ierr=0

if(is_letter(xinchar(1:1))) then   ! xinchar is a symbol
  call getparam(xinchar,cval)
else                               ! xinchar is a number string
  call value(xinchar,rval,ios)     ! rval is the value of xinchar
  if(ios > 0) then
    ierr=3
    write(*,*) 'Error: number string ',trim(xinchar),' does not correspond to a valid number' 
    write(*,*)
  end if
  cval=cmplx(rval,0.0_kr8,kc8)
  return
end if

end subroutine valuep


!**********************************************************************


SUBROUTINE PUSHOP(op)  ! Puts an operator on operator stack

type(item):: op

itop=itop+1
if(itop > numtok) then
  write(*,*) 'Error: operator stack overflow in evaluation of expression'
  write(*,*)
  return
end if
opstack(itop)=op

end subroutine pushop

SUBROUTINE POPOP(op) ! Takes top operator of operator stack and assigns it to op

type(item):: op

op=opstack(itop)
itop=itop-1

end subroutine popop

SUBROUTINE PUSHVAL(val) ! Puts value on value stack

complex(kc8) :: val

itop=itop+1
if(itop > numtok) then
  write(*,*) 'Error: value stack overflow in evaluation of expression'
  write(*,*)
  return
end if
valstack(itop)=val

end subroutine pushval

SUBROUTINE POPVAL(val) ! Takes top value off value stack and assigns it to val

complex(kc8) :: val

val=valstack(itop)
itop=itop-1

end subroutine popval

!**********************************************************************

SUBROUTINE GETPARAM_DC(sym,var)  ! Find double precision complex value var
                                 ! corresponding to symbol sym

character(len=*) :: sym
character(len=len_trim(sym)) :: usym
complex(kc8) :: var
integer ifind,j

ierr=0
sym=adjustl(sym)
if((.not.is_letter(sym(1:1))) .or. len_trim(sym)>24) then
  ierr=11
  write(*,*) 'Error: symbol ',trim(sym),' has incorrect format'
  write(*,*)
  return
end if
ifind=0
usym=uppercase(sym)
do j=1,nparams
  if(trim(usym) == trim(params(j)%symbol)) then
    var=params(j)%value
    ifind=j
    exit
  end if
end do
if(ifind == 0) then          
  ierr=4
  write(*,*) 'Error: symbol ',trim(sym), ' not in symbol table'
  write(*,*) 
  return
end if

end subroutine getparam_dc

!**********************************************************************

SUBROUTINE GETPARAM_SC(sym,var)  ! Find single precision complex value var
                                 ! corresponding to symbol sym


character(len=*) :: sym
complex(kc4) :: var
complex(kc8) :: vard

call getparam_dc(sym,vard)
var=vard

end subroutine getparam_sc

!**********************************************************************

SUBROUTINE GETPARAM_DR(sym,var)  ! Find double precision real value var
                                 ! corresponding to symbol sym


character(len=*) :: sym
real(kr8) :: var
complex(kc8) :: vard

call getparam_dc(sym,vard)
var=real(vard,kr8)

end subroutine getparam_dr

!**********************************************************************

SUBROUTINE GETPARAM_SR(sym,var)  ! Find single precision real value var
                                 ! corresponding to symbol sym


character(len=*) :: sym
real(kr4) :: var
complex(kc8) :: vard

call getparam_dc(sym,vard)
var=real(vard)

end subroutine getparam_sr

!**********************************************************************

SUBROUTINE GETPARAM_DI(sym,ivar)  ! Find double precision integer value ivar
                                  ! corresponding to symbol sym


character(len=*) :: sym
integer(ki8) :: ivar
complex(kc8) :: vard

call getparam_dc(sym,vard)
ivar=nint(real(vard,kr8),ki8)

end subroutine getparam_di

!**********************************************************************

SUBROUTINE GETPARAM_SI(sym,ivar)  ! Find single precision integer value ivar
                                  ! corresponding to symbol sym


character(len=*) :: sym
integer(ki4) :: ivar
complex(kc8) :: vard

call getparam_dc(sym,vard)
ivar=nint(real(vard,kr8),ki4)

end subroutine getparam_si

!**********************************************************************

SUBROUTINE EVALEQN(eqn)  ! Evaluate an equation

character(len=*) :: eqn
character(len=len(eqn)) :: args(2)
integer nargs

call parse(eqn,'=',args,nargs)   ! Seperate right- and left-hand-sides
call defparam(adjustl(args(1)),args(2)) ! Evaluate right-hand-side and
                                        ! assign to symbol on the
                                        ! left-hand-side.
end subroutine evaleqn

!**********************************************************************

SUBROUTINE LISTVAR      ! List all variables and their values

integer i

write(*,'(/a)') ' VARIABLE LIST:'
if(nparams == 0) then            ! Initialize symbol table
  params(1)%symbol='PI'
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
  params(2)%symbol='I'
  params(2)%value=(0.0_kr8,1.0_kr8)
  nparams=2
end if
do i=1,nparams
  write(*,*) trim(params(i)%symbol),' = ',params(i)%value
end do

end subroutine listvar

!**********************************************************************

SUBROUTINE tokenize(expr,list)    ! Evaluate expression expr for
                                    ! val double precision complex
character (len=*),intent(in) :: expr
type(Tlist) :: list

character (len=len(expr)+1) :: tempstr
integer :: isp(0:numtok)          ! On stack priority of operators in opstack
integer :: lstr
type(item):: token(numtok)      ! List of tokens ( a token is an operator or 
                                ! operand) in postfix order
type(item) :: x,junk,tok
integer i,icp,insp,isum,ntok,ifound,j
real(kr8) value

ierr=0
token(1:)%char=' '

if(nparams == 0) then                  ! Initialize symbol table
  params(1)%symbol='PI'
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
  params(2)%symbol='I'
  params(2)%value=(0.0_kr8,1.0_kr8)
  nparams=2
end if

if(len_trim(expr) == 0) then           ! Expression empty
  ierr=1
  write(*,*) 'Error: expression being evaluated is empty'
  return
end if

! ----------------------------------------------------------------------
! Check the expression does not have too many tokens
! ----------------------------------------------------------------------

     tempstr=uppercase(adjustl(expr))
     call removesp(tempstr)
     ntok=0
     do
        call get_next_token(tempstr,tok,icp,insp)
        if (tok%type.eq.'E') exit

        ntok=ntok+1
     enddo

     if (ntok.gt.numtok) then
        ierr=16
        write (*,*) 'Error: Too many tokens in expression'
        write (*,*) 'Found ',ntok,': maximum',numtok
        return
     endif

tempstr=uppercase(adjustl(expr))
call removesp(tempstr)   ! Removes spaces, tabs, and control characters

! ****************************************************************************
! STEP 1:  Convert string to token array. Each token is either an operator or
!          an operand. Token array will be in postfix (reverse Polish) order.
!*****************************************************************************

ntok=0
ibin=0
itop=0
isp(:)=0                 ! added by JFB
do
  lstr=len_trim(tempstr)
  call get_next_token(tempstr(1:lstr),tok,icp,insp)
  select case(tok%type)
  case('S')
    ntok=ntok+1
    token(ntok)=tok
  case('E')
    do 
      if(itop < 1)exit
      call popop(x)        ! Output remaining operators on stack
      ntok=ntok+1
      token(ntok)=x
    end do
    ntok=ntok+1
    token(ntok)=tok
    exit
  case('R')  ! Token is right parenenthesis
    do 
      if(opstack(itop)%type == 'L') exit  ! Output operators on stack down
      call popop(x)                       ! to left parenthesis
      ntok=ntok+1
      token(ntok)=x
    end do                             
    call popop(junk)                      ! Remove left parenthesis from stack
    if(opstack(itop)%type == 'F') then    ! Output function name if present
      call popop(x)
      ntok=ntok+1
      token(ntok)=x
    end if                               
  case('D')  ! Token is comma
    do 
      if(opstack(itop)%type == 'L') exit  ! Output operators on stack down
      call popop(x)                       ! to left parenthesis
      ntok=ntok+1
      token(ntok)=x
    end do                              
  case('U','B','L','F') ! Token is operator, left parenthesis or function name
    do 
      if(isp(itop) < icp) exit            ! Output operators on stack having
      call popop(x)                       ! an instack priority that is
      ntok=ntok+1                         ! greater than or equal to the
      token(ntok)=x                       ! priority of the incoming operator  
    end do                                
    call pushop(tok)     ! Put incoming operator on stack                       
    isp(itop)=insp
  end select
end do

isum=0                                 ! Error check for matching parentheses
do i=1,ntok
  if(token(i)%type == 'L' ) isum=isum+1
  if(token(i)%type == 'R' ) isum=isum-1
end do
if(isum /= 0) then
  ierr=2
  write(*,*) 'Error in the evaluation of the expression ',trim(expr)
  write(*,*) "Parentheses don't match"
  write(*,*)
  return
end if

list%ntok=ntok

do i=1,ntok
   x=token(i)

!   print *,i,trim(x%char)

   if (x%type.eq.'S') then
      if (is_letter(x%char(1:1))) then

         ifound=0
         do j=1,nparams
!            print *,j,trim(params(j)%symbol)
            if (params(j)%symbol.eq.x%char) then
               x%ind=j
               ifound=1
               exit
            endif
         enddo
         
         if (ifound.eq.0) then
            print *,'Tokenize_DC: could not find variable'
            return
         endif
         
         x%value=0
      else
!         print *,'x%char',trim(x%char)
         read (x%char,*) value
         x%value=value
         x%ind=0
      endif

      token(i)=x
   endif
enddo

list%token(:)=token(:)

end SUBROUTINE tokenize

! ######################################################################

subroutine evaltoken_DR(list,val)

  type(Tlist) :: list
  real(kr8) val

  complex(kr8) cval

  call evaltoken_DC(list,cval)

  val=real(cval,kind=kr8)
end subroutine evaltoken_DR

! ######################################################################

subroutine evaltoken_DC(list,val)

implicit none

type(Tlist) :: list
complex(kr8) val

character(30) :: tempstr
complex(kc8) :: cval,oper1,oper2
real(kr8) :: valr,vali

type(item) :: x
integer i,ntok

ntok=list%ntok

ierr=0
itop=0
do i=1,ntok
   x=list%token(i)

   select case(x%type)
   case('E')  ! Token is end token
      if(itop>1) then                
         ierr=12
         write(*,*) 'Error: missing operator in token list'
         write(*,*)
         return
      end if
      call popval(val)               ! Final result left on stack of values
      exit
   case('S')  ! Token is operand

      if (x%ind.eq.0) then
         cval=x%value
      else
         cval=params(x%ind)%value
      endif
     
      call pushval(cval)             ! Put value of operand on stack

   case('B')  ! Token is a binary operator
      if(itop < 2) then
         ierr=5
         write(*,*) 'Error in evaluation of token list'
         write(*,*) 'Less than two operands for binary operator  ',trim(x%char)
         write(*,*)
         return
      end if
      call popval(oper1)             ! Pull off top two values from stack
      call popval(oper2)
      select case(trim(x%char))      ! Perform operation on values
      case('^')
         cval=oper2**oper1
      case('*')
         cval=oper2*oper1
      case('/')
         if(oper1 == (0._kr8,0._kr8)) then
            ierr=10
            write(*,*) 'Error in expression (token list)'
            write(*,*) 'Division by zero'
            write(*,*)
            cval=0
         else
            cval=oper2/oper1
         end if
      case('+')
         cval=oper2+oper1
      case('-')
         cval=oper2-oper1
      end select
      call pushval(cval)             ! Put result back on stack
   case('U')  ! Token is unary operator
      if(itop == 0) then
         ierr=6
         write(*,*) 'Error in expression (token list)'
         write(*,*) 'No operand for unary operator ',trim(x%char)
         write(*,*)
         return
      else
         call popval(oper1)           ! Pull top value off stack
      end if
      select case(trim(x%char))      ! Operate on value
      case('+')
         cval=oper1
      case('-')
         cval=-oper1
      end select
      call pushval(cval)             ! Put result back on stack
   case('F')  ! Token is a function name
      if(itop == 0) then
         ierr=7
         write(*,*) 'Error in expression (token list)'
         write(*,*) 'Missing argument(s) for function ',trim(x%char)
         write(*,*)
         return
      else  
         call popval(oper1)           ! Pull top value off stack
      end if
      tempstr=uppercase(x%char)
      select case(trim(tempstr))      ! Evaluate function
      case('SIN')
         cval=sin(oper1)
      case('COS')
         cval=cos(oper1)
      case('TAN')
         oper2=cos(oper1)
         if(abs(oper2) == 0.0_kr8) then
            ierr=14
            write(*,*) 'Error: argument of tan function a multiple',&
                 ' of pi/2 in expression '
            write(*,*)
            cval=0
         else 
            cval=sin(oper1)/oper2
         endif
      case('SQRT')
         if(real(oper1,kr8) < 0. .and. aimag(oper1)==0.) then
            ierr=9
            write(*,*) 'Warning: square root of negative real number',&
                 ' in expression '
            write(*,*)
         end if
         cval=sqrt(oper1)
      case('ABS')
         cval=abs(oper1)
      case('LN')
         if(real(oper1,kr8) <= 0. .and. aimag(oper1)==0.) then
            ierr=8
            write(*,*) 'Error: negative real or zero argument for',&
                 ' natural logarithm in expression'
            write(*,*)
            cval=0
         else
            cval=log(oper1)
         end if

      case('LOG')
         if(real(oper1,kr8) <= 0. .and. aimag(oper1)==0.) then
            ierr=8
            write(*,*) 'Error: negative real or zero argument for base',&
                 '10 logarithm in expression '
            write(*,*)
            cval=0
         else
            cval=log(oper1)/2.30258509299405_kr8
         end if

      case('EXP')
         cval=exp(oper1)
      case('COMPLEX')
         if(itop == 0) then
            ierr=7
            write(*,*) 'Error in expression '
            write(*,*) 'Missing argument(s) for function ',trim(x%char)
            write(*,*)
            return
         else  
            call popval(oper2)  ! Pull second argument off stack
         end if
         valr=real(oper2,kr8)
         vali=real(oper1,kr8)
         cval=cmplx(valr,vali,kc8)
      case('CONJG')
         cval=conjg(oper1)
      case('ANG')
         cval=atan2(aimag(oper1),real(oper1,kr8))
      case('REAL')
         cval=real(oper1,kr8)
      case('IMAG')
         cval=aimag(oper1)
      case default ! Undefined function
         ierr=13
         write(*,*) 'Error: the function ',trim(x%char), ' is undefined',&
              ' in the expression '
         write(*,*)
         return
      end select
      call pushval(cval)    ! Put result back on stack
   end select
end do

end subroutine evaltoken_DC


end module evaluate