Subversion Repositories lagranto.icon

Rev

Rev 3 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 3 Rev 5
1 
module evaluate
 1 
module evaluate
2 
 
 2 
 
3 
! The user can assign values to parameters that can be used in expressions with
3 
! The user can assign values to parameters that can be used in expressions with
4 
! the subroutine defparam. The calling syntax is
 4 
! the subroutine defparam. The calling syntax is
5 
!
 5 
!
6 
!     call defparam(symbol,value) or call defparam(symbol,expr)
 6 
!     call defparam(symbol,value) or call defparam(symbol,expr)
7 
!
 7 
!
8 
! where symbol is the desired parameter name; value is a real, integer, or
8 
! where symbol is the desired parameter name; value is a real, integer, or
9 
! complex variable (single or double precision); and expr is a string
 9 
! complex variable (single or double precision); and expr is a string
10 
! containing an expression to be evaluated. The value obtained by evaluating the
 10 
! containing an expression to be evaluated. The value obtained by evaluating the
11 
! expression expr is associated with the parameter symbol. Parameter names must
 11 
! expression expr is associated with the parameter symbol. Parameter names must
12 
! begin with a letter (a-z, A-Z) and must not be longer than 24 characters.
 12 
! begin with a letter (a-z, A-Z) and must not be longer than 24 characters.
13 
! Parameter names are not case dependent.
 13 
! Parameter names are not case dependent.
14 
!
 14 
!
15 
! An expression can be evaluated with the subroutine evalexpr. The calling
 15 
! An expression can be evaluated with the subroutine evalexpr. The calling
16 
! syntax is
16 
! syntax is
17 
!
 17 
!
18 
!          call evalexpr(expr,value)
 18 
!          call evalexpr(expr,value)
19 
!
 19 
!
20 
! where expr is a string containing the expression to be evaluated; value is the
 20 
! where expr is a string containing the expression to be evaluated; value is the
21 
! result (single or double precision real, complex or integer). The
21 
! result (single or double precision real, complex or integer). The
22 
! expression can contain the arithmetic operations +, -, *, /, or ^ as well as
 22 
! expression can contain the arithmetic operations +, -, *, /, or ^ as well as
23 
! the functions sin, cos, tan, log, ln, abs, exp, sqrt, real, imag, conjg, and
23 
! the functions sin, cos, tan, log, ln, abs, exp, sqrt, real, imag, conjg, and
24 
! ang (the function ang calculates the phase angle of its complex argument). The
 24 
! ang (the function ang calculates the phase angle of its complex argument). The
25 
! expression can also contain numerical values and previously defined parameters
 25 
! expression can also contain numerical values and previously defined parameters
26 
! Grouping by nested levels of parentheses is also allowed. The parameters pi
 26 
! Grouping by nested levels of parentheses is also allowed. The parameters pi
27 
! and i (imaginary unit) are predefined. Complex numbers can be entered as a+i*b
27 
! and i (imaginary unit) are predefined. Complex numbers can be entered as a+i*b
28 
! if the parameter i has not been redefined by the user. Complex numbers can also
 28 
! if the parameter i has not been redefined by the user. Complex numbers can also
29 
! be entered using complex(a,b).
 29 
! be entered using complex(a,b).
30 
! Example expression:
30 
! Example expression:
31 
!
 31 
!
32 
!          conjg(((cos(x) + sqrt(a+i*b))^2+complex(ln(1.6e-4),20))/2)
 32 
!          conjg(((cos(x) + sqrt(a+i*b))^2+complex(ln(1.6e-4),20))/2)
33 
!
 33 
!
34 
! An equation of the form <symbol> = <expression> can be evaluated using the
34 
! An equation of the form <symbol> = <expression> can be evaluated using the
35 
! subroutine evaleqn. The calling syntax is
 35 
! subroutine evaleqn. The calling syntax is
36 
!
 36 
!
37 
!          call evaleqn(eqn)
37 
!          call evaleqn(eqn)
38 
!
 38 
!
39 
! where eqn is a string containing the equation. The right-hand-side of the
 39 
! where eqn is a string containing the equation. The right-hand-side of the
40 
! equation is evaluated and assigned to the symbol given by the left-hand-side.
 40 
! equation is evaluated and assigned to the symbol given by the left-hand-side.
41 
!
 41 
!
42 
! The value assigned to a symbol can be retrieved using the subroutine getparam.
 42 
! The value assigned to a symbol can be retrieved using the subroutine getparam.
43 
! The calling syntax is
 43 
! The calling syntax is
44 
!
 44 
!
45 
!          call getparam(sym,value)
 45 
!          call getparam(sym,value)
46 
!
 46 
!
47 
! where sym is a symbol string; value is a numeric variable (any of the six
47 
! where sym is a symbol string; value is a numeric variable (any of the six
48 
! standard types).
 48 
! standard types).
49 
!
 49 
!
50 
! The symbols and their values in the symbol table can be listed using the
 50 
! The symbols and their values in the symbol table can be listed using the
51 
! subroutine listvar. The variable ierr is always available following a call
 51 
! subroutine listvar. The variable ierr is always available following a call
52 
! to any of the above subroutines and is zero if there were no errors. The
 52 
! to any of the above subroutines and is zero if there were no errors. The
53 
! possible nonzero values for ierr are
 53 
! possible nonzero values for ierr are
54 
!
 54 
!
55 
! 1       Expression empty
 55 
! 1       Expression empty
56 
! 2       Parentheses don't match
 56 
! 2       Parentheses don't match
57 
! 3       Number string does not correspond to a valid number
 57 
! 3       Number string does not correspond to a valid number
58 
! 4       Undefined symbol
 58 
! 4       Undefined symbol
59 
! 5       Less than two operands for binary operation
 59 
! 5       Less than two operands for binary operation
60 
! 6       No operand for unary plus or minus operators
 60 
! 6       No operand for unary plus or minus operators
61 
! 7       No argument(s) for function
 61 
! 7       No argument(s) for function
62 
! 8       Zero or negative real argument for logarithm
 62 
! 8       Zero or negative real argument for logarithm
63 
! 9       Negative real argument for square root
 63 
! 9       Negative real argument for square root
64 
! 10      Division by zero
 64 
! 10      Division by zero
65 
! 11      Improper symbol format
 65 
! 11      Improper symbol format
66 
! 12      Missing operator
 66 
! 12      Missing operator
67 
! 13      Undefined function
 67 
! 13      Undefined function
68 
! 14      Argument of tangent function a multiple of pi/2
 68 
! 14      Argument of tangent function a multiple of pi/2
69 
!
 69 
!
70 
use precision
 70 
use precision
71 
use strings
 71 
use strings
72 
 
 72 
 
73 
save
 73 
save
74 
private
 74 
private
75 
public :: valuep,evalexpr,defparam,evaleqn,getparam,listvar,ierr
 75 
public :: valuep,evalexpr,defparam,evaleqn,getparam,listvar,ierr
76 
 
 76 
 
77 
type item
77 
type item
78 
  character(len=24):: char
 78 
  character(len=24):: char
79 
  character :: type
 79 
  character :: type
80 
end type item
 80 
end type item
81 
 
 81 
 
82 
type param
82 
type param
83 
  character (len=24):: symbol
 83 
  character (len=24):: symbol
84 
  complex(kc8):: value
 84 
  complex(kc8):: value
85 
end type param
 85 
end type param
86 
 
 86 
 
87 
interface defparam
87 
interface defparam
88 
  module procedure strdef       ! value given by expression
88 
  module procedure strdef       ! value given by expression
89 
  module procedure valdef_dc    ! Double precision complex value
 89 
  module procedure valdef_dc    ! Double precision complex value
90 
  module procedure valdef_sc    ! Single precision complex value
 90 
  module procedure valdef_sc    ! Single precision complex value
91 
  module procedure valdef_dr    ! Double precision real value
 91 
  module procedure valdef_dr    ! Double precision real value
92 
  module procedure valdef_sr    ! Single precision real value
 92 
  module procedure valdef_sr    ! Single precision real value
93 
  module procedure valdef_di    ! Double precision integer value
 93 
  module procedure valdef_di    ! Double precision integer value
94 
  module procedure valdef_si    ! Single precision integer value
 94 
  module procedure valdef_si    ! Single precision integer value
95 
end interface
 95 
end interface
96 
 
 96 
 
97 
interface evalexpr
 97 
interface evalexpr
98 
  module procedure evalexpr_dc  ! Double precision complex result
 98 
  module procedure evalexpr_dc  ! Double precision complex result
99 
  module procedure evalexpr_sc  ! Single precision complex result
 99 
  module procedure evalexpr_sc  ! Single precision complex result
100 
  module procedure evalexpr_dr  ! Double precision real result
 100 
  module procedure evalexpr_dr  ! Double precision real result
101 
  module procedure evalexpr_sr  ! Single precision real result
 101 
  module procedure evalexpr_sr  ! Single precision real result
102 
  module procedure evalexpr_di  ! Double precision integer result
 102 
  module procedure evalexpr_di  ! Double precision integer result
103 
  module procedure evalexpr_si  ! Single precision integer result
 103 
  module procedure evalexpr_si  ! Single precision integer result
104 
end interface
 104 
end interface
105 
 
 105 
 
106 
interface getparam
 106 
interface getparam
107 
  module procedure getparam_dc  ! Double precision complex result
 107 
  module procedure getparam_dc  ! Double precision complex result
108 
  module procedure getparam_sc  ! Single precision complex result
 108 
  module procedure getparam_sc  ! Single precision complex result
109 
  module procedure getparam_dr  ! Double precision real result
 109 
  module procedure getparam_dr  ! Double precision real result
110 
  module procedure getparam_sr  ! Single precision real result
 110 
  module procedure getparam_sr  ! Single precision real result
111 
  module procedure getparam_di  ! Double precision integer result
 111 
  module procedure getparam_di  ! Double precision integer result
112 
  module procedure getparam_si  ! Single precision integer result
 112 
  module procedure getparam_si  ! Single precision integer result
113 
end interface
 113 
end interface
114 
 
 114 
 
115 
integer,parameter :: numtok=100  ! Maximum number of tokens
115 
integer,parameter :: numtok=100  ! Maximum number of tokens
116 
type(param) :: params(100)       ! Symbol table
 116 
type(param) :: params(100)       ! Symbol table
117 
integer :: nparams=0,itop,ibin
 117 
integer :: nparams=0,itop,ibin
118 
complex(kc8) :: valstack(numtok) ! Stack used in evaluation of expression
 118 
complex(kc8) :: valstack(numtok) ! Stack used in evaluation of expression
119 
type(item):: opstack(numtok)     ! Operator stack used in conversion to postfix
 119 
type(item):: opstack(numtok)     ! Operator stack used in conversion to postfix
120 
integer :: ierr                  ! Error flag
 120 
integer :: ierr                  ! Error flag
121 
 
 121 
 
122 
!**********************************************************************
 122 
!**********************************************************************
123 
 
 123 
 
124 
contains
 124 
contains
125 
 
 125 
 
126 
!**********************************************************************
 126 
!**********************************************************************
127 
 
 127 
 
128 
 
 128 
 
129 
SUBROUTINE EVALEXPR_DC(expr,val)    ! Evaluate expression expr for
 129 
SUBROUTINE EVALEXPR_DC(expr,val)    ! Evaluate expression expr for
130 
                                    ! val double precision complex
 130 
                                    ! val double precision complex
131 
 
 131 
 
132 
character (len=*),intent(in) :: expr
 132 
character (len=*),intent(in) :: expr
133 
complex(kc8) :: val
 133 
complex(kc8) :: val
134 
character (len=len(expr)+1) :: tempstr
 134 
character (len=len(expr)+1) :: tempstr
135 
character :: cop
 135 
character :: cop
136 
integer :: isp(numtok)          ! On stack priority of operators in opstack
 136 
integer :: isp(numtok)          ! On stack priority of operators in opstack
137 
integer :: lstr
 137 
integer :: lstr
138 
complex(kc8) :: cval,oper1,oper2
 138 
complex(kc8) :: cval,oper1,oper2
139 
real(kr8) :: valr,vali
 139 
real(kr8) :: valr,vali
140 
type(item):: token(numtok)      ! List of tokens ( a token is an operator or
140 
type(item):: token(numtok)      ! List of tokens ( a token is an operator or
141 
                                ! operand) in postfix order
 141 
                                ! operand) in postfix order
142 
type(item) :: x,junk,tok
 142 
type(item) :: x,junk,tok
143 
 
 143 
 
144 
ierr=0
 144 
ierr=0
145 
token(1:)%char=' '
 145 
token(1:)%char=' '
146 
 
 146 
 
147 
if(nparams == 0) then                  ! Initialize symbol table
 147 
if(nparams == 0) then                  ! Initialize symbol table
148 
  params(1)%symbol='PI'
 148 
  params(1)%symbol='PI'
149 
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
 149 
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
150 
  params(2)%symbol='I'
 150 
  params(2)%symbol='I'
151 
  params(2)%value=(0.0_kr8,1.0_kr8)
 151 
  params(2)%value=(0.0_kr8,1.0_kr8)
152 
  nparams=2
 152 
  nparams=2
153 
end if
 153 
end if
154 
 
 154 
 
155 
if(len_trim(expr) == 0) then           ! Expression empty
 155 
if(len_trim(expr) == 0) then           ! Expression empty
156 
  ierr=1
 156 
  ierr=1
157 
  write(*,*) 'Error: expression being evaluated is empty'
 157 
  write(*,*) 'Error: expression being evaluated is empty'
158 
  return
 158 
  return
159 
end if
 159 
end if
160 
 
 160 
 
161 
tempstr=adjustl(expr)
 161 
tempstr=adjustl(expr)
162 
call removesp(tempstr)   ! Removes spaces, tabs, and control characters
 162 
call removesp(tempstr)   ! Removes spaces, tabs, and control characters
163 
 
 163 
 
164 
! ****************************************************************************
 164 
! ****************************************************************************
165 
! STEP 1:  Convert string to token array. Each token is either an operator or
 165 
! STEP 1:  Convert string to token array. Each token is either an operator or
166 
!          an operand. Token array will be in postfix (reverse Polish) order.
 166 
!          an operand. Token array will be in postfix (reverse Polish) order.
167 
!*****************************************************************************
 167 
!*****************************************************************************
168 
 
 168 
 
169 
ntok=0
 169 
ntok=0
170 
ibin=0
 170 
ibin=0
171 
itop=0
 171 
itop=0
172 
do
 172 
do
173 
  lstr=len_trim(tempstr)
 173 
  lstr=len_trim(tempstr)
174 
  call get_next_token(tempstr(1:lstr),tok,icp,insp)
 174 
  call get_next_token(tempstr(1:lstr),tok,icp,insp)
175 
  select case(tok%type)
 175 
  select case(tok%type)
176 
  case('S')
 176 
  case('S')
177 
    ntok=ntok+1
 177 
    ntok=ntok+1
178 
    token(ntok)=tok
 178 
    token(ntok)=tok
179 
  case('E')
 179 
  case('E')
180 
    do
180 
    do
181 
      if(itop < 1)exit
 181 
      if(itop < 1)exit
182 
      call popop(x)        ! Output remaining operators on stack
 182 
      call popop(x)        ! Output remaining operators on stack
183 
      ntok=ntok+1
 183 
      ntok=ntok+1
184 
      token(ntok)=x
 184 
      token(ntok)=x
185 
    end do
 185 
    end do
186 
    ntok=ntok+1
 186 
    ntok=ntok+1
187 
    token(ntok)=tok
 187 
    token(ntok)=tok
188 
    exit
 188 
    exit
189 
  case('R')  ! Token is right parenenthesis
 189 
  case('R')  ! Token is right parenenthesis
190 
    do
190 
    do
191 
      if(opstack(itop)%type == 'L') exit  ! Output operators on stack down
 191 
      if(opstack(itop)%type == 'L') exit  ! Output operators on stack down
192 
      call popop(x)                       ! to left parenthesis
 192 
      call popop(x)                       ! to left parenthesis
193 
      ntok=ntok+1
 193 
      ntok=ntok+1
194 
      token(ntok)=x
 194 
      token(ntok)=x
195 
    end do
195 
    end do
196 
    call popop(junk)                      ! Remove left parenthesis from stack
 196 
    call popop(junk)                      ! Remove left parenthesis from stack
197 
    if(opstack(itop)%type == 'F') then    ! Output function name if present
 197 
    if(opstack(itop)%type == 'F') then    ! Output function name if present
198 
      call popop(x)
 198 
      call popop(x)
199 
      ntok=ntok+1
 199 
      ntok=ntok+1
200 
      token(ntok)=x
 200 
      token(ntok)=x
201 
    end if
201 
    end if
202 
  case('D')  ! Token is comma
 202 
  case('D')  ! Token is comma
203 
    do
203 
    do
204 
      if(opstack(itop)%type == 'L') exit  ! Output operators on stack down
 204 
      if(opstack(itop)%type == 'L') exit  ! Output operators on stack down
205 
      call popop(x)                       ! to left parenthesis
 205 
      call popop(x)                       ! to left parenthesis
206 
      ntok=ntok+1
 206 
      ntok=ntok+1
207 
      token(ntok)=x
 207 
      token(ntok)=x
208 
    end do
208 
    end do
209 
  case('U','B','L','F') ! Token is operator, left parenthesis or function name
 209 
  case('U','B','L','F') ! Token is operator, left parenthesis or function name
210 
    do
210 
    do
211 
      if(isp(itop) < icp) exit            ! Output operators on stack having
 211 
      if(isp(itop) < icp) exit            ! Output operators on stack having
212 
      call popop(x)                       ! an instack priority that is
 212 
      call popop(x)                       ! an instack priority that is
213 
      ntok=ntok+1                         ! greater than or equal to the
 213 
      ntok=ntok+1                         ! greater than or equal to the
214 
      token(ntok)=x                       ! priority of the incoming operator
214 
      token(ntok)=x                       ! priority of the incoming operator
215 
    end do
215 
    end do
216 
    call pushop(tok)     ! Put incoming operator on stack
216 
    call pushop(tok)     ! Put incoming operator on stack
217 
    isp(itop)=insp
 217 
    isp(itop)=insp
218 
  end select
 218 
  end select
219 
end do
 219 
end do
220 
 
 220 
 
221 
isum=0                                 ! Error check for matching parentheses
 221 
isum=0                                 ! Error check for matching parentheses
222 
do i=1,ntok
 222 
do i=1,ntok
223 
  if(token(i)%type == 'L' ) isum=isum+1
 223 
  if(token(i)%type == 'L' ) isum=isum+1
224 
  if(token(i)%type == 'R' ) isum=isum-1
 224 
  if(token(i)%type == 'R' ) isum=isum-1
225 
end do
 225 
end do
226 
if(isum /= 0) then
 226 
if(isum /= 0) then
227 
  ierr=2
 227 
  ierr=2
228 
  write(*,*) 'Error in the evaluation of the expression ',trim(expr)
 228 
  write(*,*) 'Error in the evaluation of the expression ',trim(expr)
229 
  write(*,*) "Parentheses don't match"
 229 
  write(*,*) "Parentheses don't match"
230 
  write(*,*)
 230 
  write(*,*)
231 
  return
 231 
  return
232 
end if
 232 
end if
233 
 
 233 
 
234 
 
 234 
 
235 
!*****************************************************************************
 235 
!*****************************************************************************
236 
! STEP 2: Evaluate token string in postfix order
 236 
! STEP 2: Evaluate token string in postfix order
237 
!*****************************************************************************
 237 
!*****************************************************************************
238 
 
 238 
 
239 
itop=0
 239 
itop=0
240 
do i=1,ntok
 240 
do i=1,ntok
241 
  x=token(i)
 241 
  x=token(i)
242 
  select case(x%type)
 242 
  select case(x%type)
243 
  case('E')  ! Token is end token
 243 
  case('E')  ! Token is end token
244 
    if(itop>1) then
244 
    if(itop>1) then
245 
      ierr=12
 245 
      ierr=12
246 
      write(*,*) 'Error: missing operator in expression ',trim(expr)
 246 
      write(*,*) 'Error: missing operator in expression ',trim(expr)
247 
      write(*,*)
 247 
      write(*,*)
248 
      return
 248 
      return
249 
    end if
 249 
    end if
250 
    call popval(val)               ! Final result left on stack of values
 250 
    call popval(val)               ! Final result left on stack of values
251 
    exit
 251 
    exit
252 
  case('S')  ! Token is operand
 252 
  case('S')  ! Token is operand
253 
    call valuep(x%char,cval)       ! Evaluate operand
 253 
    call valuep(x%char,cval)       ! Evaluate operand
254 
    if(ierr/=0) return
 254 
    if(ierr/=0) return
255 
    call pushval(cval)             ! Put value of operand on stack
 255 
    call pushval(cval)             ! Put value of operand on stack
256 
  case('B')  ! Token is a binary operator
 256 
  case('B')  ! Token is a binary operator
257 
    if(itop < 2) then
 257 
    if(itop < 2) then
258 
      ierr=5
 258 
      ierr=5
259 
      write(*,*) 'Error in evaluation of expression ',trim(expr)
 259 
      write(*,*) 'Error in evaluation of expression ',trim(expr)
260 
      write(*,*) 'Less than two operands for binary operator  '&
 260 
      write(*,*) 'Less than two operands for binary operator  '&
261 
                 ,trim(x%char)
 261 
                 ,trim(x%char)
262 
      write(*,*)
 262 
      write(*,*)
263 
      return
 263 
      return
264 
    end if
264 
    end if
265 
    call popval(oper1)             ! Pull off top two values from stack
 265 
    call popval(oper1)             ! Pull off top two values from stack
266 
    call popval(oper2)
 266 
    call popval(oper2)
267 
    select case(trim(x%char))      ! Perform operation on values
 267 
    select case(trim(x%char))      ! Perform operation on values
268 
    case('^')
 268 
    case('^')
269 
      cval=oper2**oper1
 269 
      cval=oper2**oper1
270 
    case('*')
 270 
    case('*')
271 
      cval=oper2*oper1
 271 
      cval=oper2*oper1
272 
    case('/')
 272 
    case('/')
273 
      if(oper1 == (0._kr8,0._kr8)) then
 273 
      if(oper1 == (0._kr8,0._kr8)) then
274 
        ierr=10
 274 
        ierr=10
275 
        write(*,*) 'Error in expression ',trim(expr)
 275 
        write(*,*) 'Error in expression ',trim(expr)
276 
        write(*,*) 'Division by zero'
 276 
        write(*,*) 'Division by zero'
277 
        write(*,*)
 277 
        write(*,*)
278 
        return
 278 
        return
279 
      end if
 279 
      end if
280 
      cval=oper2/oper1
 280 
      cval=oper2/oper1
281 
    case('+')
 281 
    case('+')
282 
      cval=oper2+oper1
 282 
      cval=oper2+oper1
283 
    case('-')
 283 
    case('-')
284 
      cval=oper2-oper1
 284 
      cval=oper2-oper1
285 
    end select
 285 
    end select
286 
    call pushval(cval)             ! Put result back on stack
 286 
    call pushval(cval)             ! Put result back on stack
287 
  case('U')  ! Token is unary operator
 287 
  case('U')  ! Token is unary operator
288 
    if(itop == 0) then
 288 
    if(itop == 0) then
289 
      ierr=6
 289 
      ierr=6
290 
      write(*,*) 'Error in expression ',trim(expr)
 290 
      write(*,*) 'Error in expression ',trim(expr)
291 
      write(*,*) 'No operand for unary operator ',trim(x%char)
 291 
      write(*,*) 'No operand for unary operator ',trim(x%char)
292 
      write(*,*)
 292 
      write(*,*)
293 
      return
 293 
      return
294 
    else
 294 
    else
295 
      call popval(oper1)           ! Pull top value off stack
 295 
      call popval(oper1)           ! Pull top value off stack
296 
    end if
 296 
    end if
297 
    select case(trim(x%char))      ! Operate on value
 297 
    select case(trim(x%char))      ! Operate on value
298 
    case('+')
 298 
    case('+')
299 
      cval=oper1
 299 
      cval=oper1
300 
    case('-')
 300 
    case('-')
301 
      cval=-oper1
 301 
      cval=-oper1
302 
    end select
 302 
    end select
303 
    call pushval(cval)             ! Put result back on stack
 303 
    call pushval(cval)             ! Put result back on stack
304 
  case('F')  ! Token is a function name
 304 
  case('F')  ! Token is a function name
305 
    if(itop == 0) then
 305 
    if(itop == 0) then
306 
      ierr=7
 306 
      ierr=7
307 
      write(*,*) 'Error in expression ',trim(expr)
 307 
      write(*,*) 'Error in expression ',trim(expr)
308 
      write(*,*) 'Missing argument(s) for function ',trim(x%char)
 308 
      write(*,*) 'Missing argument(s) for function ',trim(x%char)
309 
      write(*,*)
 309 
      write(*,*)
310 
      return
 310 
      return
311 
    else
311 
    else
312 
      call popval(oper1)           ! Pull top value off stack
 312 
      call popval(oper1)           ! Pull top value off stack
313 
    end if
313 
    end if
314 
    tempstr=uppercase(x%char)
 314 
    tempstr=uppercase(x%char)
315 
    select case(trim(tempstr))      ! Evaluate function
 315 
    select case(trim(tempstr))      ! Evaluate function
316 
    case('SIN')
 316 
    case('SIN')
317 
      cval=sin(oper1)
 317 
      cval=sin(oper1)
318 
    case('COS')
 318 
    case('COS')
319 
      cval=cos(oper1)
 319 
      cval=cos(oper1)
320 
    case('TAN')
 320 
    case('TAN')
321 
      oper2=cos(oper1)
 321 
      oper2=cos(oper1)
322 
      if(abs(oper2) == 0.0_kr8) then
 322 
      if(abs(oper2) == 0.0_kr8) then
323 
        ierr=14
 323 
        ierr=14
324 
        write(*,*) 'Error: argument of tan function a multiple',&
 324 
        write(*,*) 'Error: argument of tan function a multiple',&
325 
        ' of pi/2 in expression ',trim(expr)
 325 
        ' of pi/2 in expression ',trim(expr)
326 
        write(*,*)
 326 
        write(*,*)
327 
        return
 327 
        return
328 
      else
328 
      else
329 
        cval=sin(oper1)/oper2
 329 
        cval=sin(oper1)/oper2
330 
      endif
 330 
      endif
331 
    case('SQRT')
 331 
    case('SQRT')
332 
      if(real(oper1,kr8) < 0. .and. aimag(oper1)==0.) then
 332 
      if(real(oper1,kr8) < 0. .and. aimag(oper1)==0.) then
333 
        ierr=9
 333 
        ierr=9
334 
        write(*,*) 'Warning: square root of negative real number',&
 334 
        write(*,*) 'Warning: square root of negative real number',&
335 
                   ' in expression ',trim(expr)
 335 
                   ' in expression ',trim(expr)
336 
        write(*,*)
 336 
        write(*,*)
337 
      end if
 337 
      end if
338 
      cval=sqrt(oper1)
 338 
      cval=sqrt(oper1)
339 
    case('ABS')
 339 
    case('ABS')
340 
      cval=abs(oper1)
 340 
      cval=abs(oper1)
341 
    case('LN')
 341 
    case('LN')
342 
      if(real(oper1,kr8) <= 0. .and. aimag(oper1)==0.) then
 342 
      if(real(oper1,kr8) <= 0. .and. aimag(oper1)==0.) then
343 
        ierr=8
 343 
        ierr=8
344 
        write(*,*) 'Error: negative real or zero argument for',&
 344 
        write(*,*) 'Error: negative real or zero argument for',&
345 
                   ' natural logarithm in expression ',trim(expr)
 345 
                   ' natural logarithm in expression ',trim(expr)
346 
        write(*,*)
 346 
        write(*,*)
347 
        return
 347 
        return
348 
      end if
 348 
      end if
349 
      cval=log(oper1)
 349 
      cval=log(oper1)
350 
    case('LOG')
 350 
    case('LOG')
351 
      if(real(oper1,kr8) <= 0. .and. aimag(oper1)==0.) then
 351 
      if(real(oper1,kr8) <= 0. .and. aimag(oper1)==0.) then
352 
        ierr=8
 352 
        ierr=8
353 
        write(*,*) 'Error: negative real or zero argument for base',&
 353 
        write(*,*) 'Error: negative real or zero argument for base',&
354 
                   '10 logarithm in expression ',trim(expr)
 354 
                   '10 logarithm in expression ',trim(expr)
355 
        write(*,*)
 355 
        write(*,*)
356 
        return
 356 
        return
357 
      end if
 357 
      end if
358 
      cval=log(oper1)/2.30258509299405_kr8
 358 
      cval=log(oper1)/2.30258509299405_kr8
359 
    case('EXP')
 359 
    case('EXP')
360 
      cval=exp(oper1)
 360 
      cval=exp(oper1)
361 
    case('COMPLEX')
 361 
    case('COMPLEX')
362 
      if(itop == 0) then
 362 
      if(itop == 0) then
363 
        ierr=7
 363 
        ierr=7
364 
        write(*,*) 'Error in expression ',trim(expr)
 364 
        write(*,*) 'Error in expression ',trim(expr)
365 
        write(*,*) 'Missing argument(s) for function ',trim(x%char)
 365 
        write(*,*) 'Missing argument(s) for function ',trim(x%char)
366 
        write(*,*)
 366 
        write(*,*)
367 
        return
 367 
        return
368 
      else
368 
      else
369 
        call popval(oper2)  ! Pull second argument off stack
 369 
        call popval(oper2)  ! Pull second argument off stack
370 
      end if
370 
      end if
371 
      valr=real(oper2,kr8)
 371 
      valr=real(oper2,kr8)
372 
      vali=real(oper1,kr8)
 372 
      vali=real(oper1,kr8)
373 
      cval=cmplx(valr,vali,kc8)
 373 
      cval=cmplx(valr,vali,kc8)
374 
    case('CONJG')
 374 
    case('CONJG')
375 
      cval=conjg(oper1)
 375 
      cval=conjg(oper1)
376 
    case('ANG')
 376 
    case('ANG')
377 
      cval=atan2(aimag(oper1),real(oper1,kr8))
 377 
      cval=atan2(aimag(oper1),real(oper1,kr8))
378 
    case('REAL')
 378 
    case('REAL')
379 
      cval=real(oper1,kr8)
 379 
      cval=real(oper1,kr8)
380 
    case('IMAG')
 380 
    case('IMAG')
381 
      cval=aimag(oper1)
 381 
      cval=aimag(oper1)
382 
    case default ! Undefined function
 382 
    case default ! Undefined function
383 
      ierr=13
 383 
      ierr=13
384 
      write(*,*) 'Error: the function ',trim(x%char), ' is undefined',&
 384 
      write(*,*) 'Error: the function ',trim(x%char), ' is undefined',&
385 
                 ' in the expression ',trim(expr)
 385 
                 ' in the expression ',trim(expr)
386 
      write(*,*)
 386 
      write(*,*)
387 
      return
 387 
      return
388 
    end select
 388 
    end select
389 
    call pushval(cval)    ! Put result back on stack
 389 
    call pushval(cval)    ! Put result back on stack
390 
  end select
 390 
  end select
391 
end do
 391 
end do
392 
 
 392 
 
393 
end subroutine evalexpr_dc
 393 
end subroutine evalexpr_dc
394 
 
 394 
 
395 
!**********************************************************************
 395 
!**********************************************************************
396 
 
 396 
 
397 
SUBROUTINE GET_NEXT_TOKEN(str,tok,icp,isp)
 397 
SUBROUTINE GET_NEXT_TOKEN(str,tok,icp,isp)
398 
 
 398 
 
399 
character(len=*) :: str
 399 
character(len=*) :: str
400 
character :: cop,chtemp
 400 
character :: cop,chtemp
401 
type(item) :: tok
 401 
type(item) :: tok
402 
integer :: icp
 402 
integer :: icp
403 
 
 403 
 
404 
lstr=len_trim(str)
 404 
lstr=len_trim(str)
405 
if(lstr == 0) then
 405 
if(lstr == 0) then
406 
  tok%char='#'             ! Output end token
 406 
  tok%char='#'             ! Output end token
407 
  tok%type='E'
 407 
  tok%type='E'
408 
  return
 408 
  return
409 
end if
 409 
end if
410 
ipos=scan(str,'+-*/^(),')  ! Look for an arithmetic operator
410 
ipos=scan(str,'+-*/^(),')  ! Look for an arithmetic operator
411 
                           ! + - * / ^ ( ) or ,
 411 
                           ! + - * / ^ ( ) or ,
412 
cop=str(ipos:ipos)
412 
cop=str(ipos:ipos)
413 
select case (ipos)
413 
select case (ipos)
414 
case(0)    ! Operators not present
 414 
case(0)    ! Operators not present
415 
  ntok=ntok+1
 415 
  ntok=ntok+1
416 
  tok%char=str
 416 
  tok%char=str
417 
  tok%type='S'
 417 
  tok%type='S'
418 
  str=''
 418 
  str=''
419 
  icp=0
 419 
  icp=0
420 
  isp=0
 420 
  isp=0
421 
case(1)
421 
case(1)
422 
  tok%char=cop
 422 
  tok%char=cop
423 
  select case(cop)
 423 
  select case(cop)
424 
  case('+','-')
 424 
  case('+','-')
425 
    if(ibin==0) then
 425 
    if(ibin==0) then
426 
      tok%type='U'
 426 
      tok%type='U'
427 
      icp=4
 427 
      icp=4
428 
      isp=3
 428 
      isp=3
429 
    else
 429 
    else
430 
      tok%type='B'
 430 
      tok%type='B'
431 
      icp=1
 431 
      icp=1
432 
      isp=1
 432 
      isp=1
433 
    end if
 433 
    end if
434 
    ibin=0
 434 
    ibin=0
435 
  case('*','/')
 435 
  case('*','/')
436 
    tok%type='B'
 436 
    tok%type='B'
437 
    icp=2
 437 
    icp=2
438 
    isp=2
 438 
    isp=2
439 
    ibin=0
 439 
    ibin=0
440 
  case('^')
 440 
  case('^')
441 
    tok%type='B'
 441 
    tok%type='B'
442 
    icp=4
 442 
    icp=4
443 
    isp=3
 443 
    isp=3
444 
    ibin=0
 444 
    ibin=0
445 
  case('(')
 445 
  case('(')
446 
    tok%type='L'
 446 
    tok%type='L'
447 
    icp=4
 447 
    icp=4
448 
    isp=0
 448 
    isp=0
449 
    ibin=0
 449 
    ibin=0
450 
  case(')')
 450 
  case(')')
451 
    tok%type='R'
 451 
    tok%type='R'
452 
    icp=0
 452 
    icp=0
453 
    isp=0
 453 
    isp=0
454 
    ibin=1
 454 
    ibin=1
455 
  case(',')
 455 
  case(',')
456 
    tok%type='D'
 456 
    tok%type='D'
457 
    icp=0
 457 
    icp=0
458 
    isp=0
 458 
    isp=0
459 
    ibin=0
 459 
    ibin=0
460 
  end select
 460 
  end select
461 
  str=str(2:)
 461 
  str=str(2:)
462 
case(2:)
 462 
case(2:)
463 
  select case(cop)
 463 
  select case(cop)
464 
  case('(')
 464 
  case('(')
465 
    tok%char=str(1:ipos-1)
 465 
    tok%char=str(1:ipos-1)
466 
    tok%type='F'
 466 
    tok%type='F'
467 
    icp=4
 467 
    icp=4
468 
    isp=0
 468 
    isp=0
469 
    ibin=0
 469 
    ibin=0
470 
    str=str(ipos:)
 470 
    str=str(ipos:)
471 
  case('+','-')
 471 
  case('+','-')
472 
    chtemp=uppercase(str(ipos-1:ipos-1))
 472 
    chtemp=uppercase(str(ipos-1:ipos-1))
473 
    if(is_letter(str(1:1))==.true. .or. chtemp/='E') then
 473 
    if(is_letter(str(1:1)).eqv..true. .or. chtemp/='E') then
474 
      tok%char=str(1:ipos-1)
 474 
      tok%char=str(1:ipos-1)
475 
      tok%type='S'
 475 
      tok%type='S'
476 
      icp=0
 476 
      icp=0
477 
      isp=0
 477 
      isp=0
478 
      ibin=1
 478 
      ibin=1
479 
      str=str(ipos:)
 479 
      str=str(ipos:)
480 
    else
 480 
    else
481 
      inext=scan(str(ipos+1:),'+-*/^(),')
 481 
      inext=scan(str(ipos+1:),'+-*/^(),')
482 
      if(inext==0) then
 482 
      if(inext==0) then
483 
        tok%char=str
 483 
        tok%char=str
484 
        tok%type='S'
 484 
        tok%type='S'
485 
        icp=0
 485 
        icp=0
486 
        isp=0
 486 
        isp=0
487 
        ibin=0
 487 
        ibin=0
488 
        str=''
 488 
        str=''
489 
      else
 489 
      else
490 
        tok%char=str(1:ipos+inext-1)
 490 
        tok%char=str(1:ipos+inext-1)
491 
        tok%type='S'
 491 
        tok%type='S'
492 
        icp=0
 492 
        icp=0
493 
        isp=0
 493 
        isp=0
494 
        ibin=1
 494 
        ibin=1
495 
        str=str(ipos+inext:)
 495 
        str=str(ipos+inext:)
496 
      end if
 496 
      end if
497 
    end if
 497 
    end if
498 
  case default
 498 
  case default
499 
    tok%char=str(1:ipos-1)
 499 
    tok%char=str(1:ipos-1)
500 
    tok%type='S'
 500 
    tok%type='S'
501 
    icp=0
 501 
    icp=0
502 
    isp=0
 502 
    isp=0
503 
    ibin=1
 503 
    ibin=1
504 
    str=str(ipos:)
 504 
    str=str(ipos:)
505 
  end select
 505 
  end select
506 
end select
 506 
end select
507 
 
 507 
 
508 
end subroutine get_next_token
 508 
end subroutine get_next_token
509 
 
 509 
 
510 
 
 510 
 
511 
!**********************************************************************
 511 
!**********************************************************************
512 
 
 512 
 
513 
SUBROUTINE EVALEXPR_SC(expr,val)    ! Evaluate expression expr for
 513 
SUBROUTINE EVALEXPR_SC(expr,val)    ! Evaluate expression expr for
514 
                                    ! val single precision complex
 514 
                                    ! val single precision complex
515 
character(len=*) :: expr
 515 
character(len=*) :: expr
516 
complex(kc4) :: val
 516 
complex(kc4) :: val
517 
complex(kc8) :: vald
 517 
complex(kc8) :: vald
518 
 
 518 
 
519 
call evalexpr_dc(expr,vald)
 519 
call evalexpr_dc(expr,vald)
520 
val=vald
 520 
val=vald
521 
 
 521 
 
522 
end subroutine evalexpr_sc
 522 
end subroutine evalexpr_sc
523 
 
 523 
 
524 
!**********************************************************************
 524 
!**********************************************************************
525 
 
 525 
 
526 
SUBROUTINE EVALEXPR_SR(expr,val)    ! Evaluate expression expr for
 526 
SUBROUTINE EVALEXPR_SR(expr,val)    ! Evaluate expression expr for
527 
                                    ! val single precision real
 527 
                                    ! val single precision real
528 
character(len=*) :: expr
 528 
character(len=*) :: expr
529 
real(kr4) :: val
 529 
real(kr4) :: val
530 
complex(kc8) :: vald
 530 
complex(kc8) :: vald
531 
 
 531 
 
532 
call evalexpr_dc(expr,vald)
 532 
call evalexpr_dc(expr,vald)
533 
val=real(vald)
 533 
val=real(vald)
534 
 
 534 
 
535 
end subroutine evalexpr_sr
 535 
end subroutine evalexpr_sr
536 
 
 536 
 
537 
!**********************************************************************
 537 
!**********************************************************************
538 
 
 538 
 
539 
SUBROUTINE EVALEXPR_DR(expr,val)    ! Evaluate expression expr for
539 
SUBROUTINE EVALEXPR_DR(expr,val)    ! Evaluate expression expr for
540 
                                    ! val double precision real
 540 
                                    ! val double precision real
541 
character(len=*) :: expr
 541 
character(len=*) :: expr
542 
real(kr8) :: val
 542 
real(kr8) :: val
543 
complex(kc8) :: vald
 543 
complex(kc8) :: vald
544 
 
 544 
 
545 
call evalexpr_dc(expr,vald)
 545 
call evalexpr_dc(expr,vald)
546 
val=real(vald,kr8)
 546 
val=real(vald,kr8)
547 
 
 547 
 
548 
end subroutine evalexpr_dr
 548 
end subroutine evalexpr_dr
549 
 
 549 
 
550 
!**********************************************************************
 550 
!**********************************************************************
551 
 
 551 
 
552 
SUBROUTINE EVALEXPR_SI(expr,ival)   ! Evaluate expression expr for
552 
SUBROUTINE EVALEXPR_SI(expr,ival)   ! Evaluate expression expr for
553 
                                    ! ival single precision integer
 553 
                                    ! ival single precision integer
554 
character(len=*) :: expr
 554 
character(len=*) :: expr
555 
integer(ki4) :: ival
 555 
integer(ki4) :: ival
556 
complex(kc8) :: vald
 556 
complex(kc8) :: vald
557 
 
 557 
 
558 
call evalexpr_dc(expr,vald)
 558 
call evalexpr_dc(expr,vald)
559 
ival=nint(real(vald,kr8),ki4)
 559 
ival=nint(real(vald,kr8),ki4)
560 
 
 560 
 
561 
end subroutine evalexpr_si
 561 
end subroutine evalexpr_si
562 
 
 562 
 
563 
!**********************************************************************
 563 
!**********************************************************************
564 
 
 564 
 
565 
SUBROUTINE EVALEXPR_DI(expr,ival)   ! Evaluate expression expr for
565 
SUBROUTINE EVALEXPR_DI(expr,ival)   ! Evaluate expression expr for
566 
                                    ! ival double precision integer
 566 
                                    ! ival double precision integer
567 
character(len=*) :: expr
 567 
character(len=*) :: expr
568 
integer(ki8) :: ival
 568 
integer(ki8) :: ival
569 
complex(kc8) :: vald
 569 
complex(kc8) :: vald
570 
 
 570 
 
571 
call evalexpr_dc(expr,vald)
 571 
call evalexpr_dc(expr,vald)
572 
ival=nint(real(vald,kr8),ki8)
 572 
ival=nint(real(vald,kr8),ki8)
573 
 
 573 
 
574 
end subroutine evalexpr_di
 574 
end subroutine evalexpr_di
575 
 
 575 
 
576 
 
 576 
 
577 
!**********************************************************************
 577 
!**********************************************************************
578 
SUBROUTINE VALDEF_DC(sym,val)    ! Associates sym with val in symbol table,
 578 
SUBROUTINE VALDEF_DC(sym,val)    ! Associates sym with val in symbol table,
579 
                                 ! val double precision complex
 579 
                                 ! val double precision complex
580 
character(len=*) :: sym
 580 
character(len=*) :: sym
581 
character(len=len_trim(sym)) :: usym
 581 
character(len=len_trim(sym)) :: usym
582 
complex(kc8) :: val
 582 
complex(kc8) :: val
583 
 
 583 
 
584 
ierr=0
 584 
ierr=0
585 
if(nparams == 0) then               ! Initialize symbol table
 585 
if(nparams == 0) then               ! Initialize symbol table
586 
  params(1)%symbol='PI'
 586 
  params(1)%symbol='PI'
587 
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
 587 
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
588 
  params(2)%symbol='I'
 588 
  params(2)%symbol='I'
589 
  params(2)%value=(0.0_kr8,1.0_kr8)
 589 
  params(2)%value=(0.0_kr8,1.0_kr8)
590 
  nparams=2
 590 
  nparams=2
591 
end if
 591 
end if
592 
 
 592 
 
593 
! Assign val to sym if sym is already in symbol table
 593 
! Assign val to sym if sym is already in symbol table
594 
usym=uppercase(sym)
 594 
usym=uppercase(sym)
595 
if(is_letter(sym(1:1))==.false. .or. len_trim(sym)>24) then
 595 
if(is_letter(sym(1:1)).eqv..false. .or. len_trim(sym)>24) then
596 
  ierr=11
 596 
  ierr=11
597 
  write(*,*) 'Error: symbol ',trim(sym),' has improper format'
 597 
  write(*,*) 'Error: symbol ',trim(sym),' has improper format'
598 
  write(*,*)
 598 
  write(*,*)
599 
  return
 599 
  return
600 
end if
 600 
end if
601 
do i=1,nparams
 601 
do i=1,nparams
602 
  if(trim(usym)==trim(params(i)%symbol)) then
 602 
  if(trim(usym)==trim(params(i)%symbol)) then
603 
    params(i)%value=val
 603 
    params(i)%value=val
604 
    return
 604 
    return
605 
  end if
 605 
  end if
606 
end do
 606 
end do
607 
 
 607 
 
608 
nparams=nparams+1    ! Otherwise assign val to new symbol sym
 608 
nparams=nparams+1    ! Otherwise assign val to new symbol sym
609 
params(nparams)%symbol=usym
 609 
params(nparams)%symbol=usym
610 
params(nparams)%value=val
 610 
params(nparams)%value=val
611 
 
 611 
 
612 
end subroutine valdef_dc
 612 
end subroutine valdef_dc
613 
 
 613 
 
614 
 
 614 
 
615 
!**********************************************************************
 615 
!**********************************************************************
616 
 
 616 
 
617 
SUBROUTINE VALDEF_SC(sym,val)     ! Associates sym with val in symbol table,
 617 
SUBROUTINE VALDEF_SC(sym,val)     ! Associates sym with val in symbol table,
618 
                                  ! val single precision complex
 618 
                                  ! val single precision complex
619 
character(len=*) :: sym
 619 
character(len=*) :: sym
620 
complex(kc4) :: val
 620 
complex(kc4) :: val
621 
complex(kc8) :: vald
 621 
complex(kc8) :: vald
622 
 
 622 
 
623 
vald=val
 623 
vald=val
624 
call valdef_dc(sym,vald)
 624 
call valdef_dc(sym,vald)
625 
 
 625 
 
626 
end subroutine valdef_sc
 626 
end subroutine valdef_sc
627 
 
 627 
 
628 
 
 628 
 
629 
!**********************************************************************
 629 
!**********************************************************************
630 
 
 630 
 
631 
SUBROUTINE VALDEF_DR(sym,val)    ! Associates sym with val in symbol table,
 631 
SUBROUTINE VALDEF_DR(sym,val)    ! Associates sym with val in symbol table,
632 
                                 ! val double precision real
 632 
                                 ! val double precision real
633 
character(len=*) :: sym
 633 
character(len=*) :: sym
634 
real(kr8) :: val
 634 
real(kr8) :: val
635 
complex(kc8) :: vald
 635 
complex(kc8) :: vald
636 
 
 636 
 
637 
vald=cmplx(val,0.0_kr8,kc8)
 637 
vald=cmplx(val,0.0_kr8,kc8)
638 
call valdef_dc(sym,vald)
 638 
call valdef_dc(sym,vald)
639 
 
 639 
 
640 
end subroutine valdef_dr
 640 
end subroutine valdef_dr
641 
 
 641 
 
642 
 
 642 
 
643 
!**********************************************************************
 643 
!**********************************************************************
644 
 
 644 
 
645 
SUBROUTINE VALDEF_SR(sym,val)    ! Associates sym with val in symbol table,
 645 
SUBROUTINE VALDEF_SR(sym,val)    ! Associates sym with val in symbol table,
646 
                                 ! val single precision real
 646 
                                 ! val single precision real
647 
character(len=*) :: sym
 647 
character(len=*) :: sym
648 
real(kr4) :: val
 648 
real(kr4) :: val
649 
complex(kc8) :: vald
 649 
complex(kc8) :: vald
650 
 
 650 
 
651 
vald=cmplx(val,0.0,kc8)
 651 
vald=cmplx(val,0.0,kc8)
652 
call valdef_dc(sym,vald)
 652 
call valdef_dc(sym,vald)
653 
 
 653 
 
654 
end subroutine valdef_sr
 654 
end subroutine valdef_sr
655 
 
 655 
 
656 
 
 656 
 
657 
!**********************************************************************
 657 
!**********************************************************************
658 
 
 658 
 
659 
SUBROUTINE VALDEF_DI(sym,ival)   ! Associates sym with ival in symbol table,
 659 
SUBROUTINE VALDEF_DI(sym,ival)   ! Associates sym with ival in symbol table,
660 
                                 ! ival double precision integer
660 
                                 ! ival double precision integer
661 
character(len=*) :: sym
 661 
character(len=*) :: sym
662 
integer(ki8) :: ival
 662 
integer(ki8) :: ival
663 
complex(kc8) :: vald
 663 
complex(kc8) :: vald
664 
 
 664 
 
665 
vald=cmplx(real(ival,kr8),0.0_kr8,kc8)
 665 
vald=cmplx(real(ival,kr8),0.0_kr8,kc8)
666 
call valdef_dc(sym,vald)
 666 
call valdef_dc(sym,vald)
667 
 
 667 
 
668 
end subroutine valdef_di
 668 
end subroutine valdef_di
669 
 
 669 
 
670 
 
 670 
 
671 
!**********************************************************************
 671 
!**********************************************************************
672 
 
 672 
 
673 
SUBROUTINE VALDEF_SI(sym,ival)   ! Associates sym with ival in symbol table,
 673 
SUBROUTINE VALDEF_SI(sym,ival)   ! Associates sym with ival in symbol table,
674 
                                 ! ival single precision integer
 674 
                                 ! ival single precision integer
675 
character(len=*) :: sym
 675 
character(len=*) :: sym
676 
integer(ki4) :: ival
 676 
integer(ki4) :: ival
677 
complex(kc8) :: vald
 677 
complex(kc8) :: vald
678 
 
 678 
 
679 
vald=cmplx(real(ival,kr8),0.0,kc8)
 679 
vald=cmplx(real(ival,kr8),0.0,kc8)
680 
call valdef_dc(sym,vald)
 680 
call valdef_dc(sym,vald)
681 
 
 681 
 
682 
end subroutine valdef_si
 682 
end subroutine valdef_si
683 
 
 683 
 
684 
 
 684 
 
685 
!**********************************************************************
 685 
!**********************************************************************
686 
 
 686 
 
687 
SUBROUTINE STRDEF(sym,expr)      ! Associates sym with the value of the
 687 
SUBROUTINE STRDEF(sym,expr)      ! Associates sym with the value of the
688 
                                 ! expression expr
 688 
                                 ! expression expr
689 
 
 689 
 
690 
character(len=*) :: sym,expr
 690 
character(len=*) :: sym,expr
691 
complex(kc8) :: val
 691 
complex(kc8) :: val
692 
 
 692 
 
693 
if(nparams == 0) then            ! Initialize symbol table
 693 
if(nparams == 0) then            ! Initialize symbol table
694 
  params(1)%symbol='PI'
 694 
  params(1)%symbol='PI'
695 
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
 695 
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
696 
  params(2)%symbol='I'
 696 
  params(2)%symbol='I'
697 
  params(2)%value=(0.0_kr8,1.0_kr8)
 697 
  params(2)%value=(0.0_kr8,1.0_kr8)
698 
  nparams=2
 698 
  nparams=2
699 
end if
 699 
end if
700 
 
 700 
 
701 
call evalexpr_dc(expr,val)       ! val is value of expression expr
 701 
call evalexpr_dc(expr,val)       ! val is value of expression expr
702 
if(ierr==0 .or. ierr==9) then
 702 
if(ierr==0 .or. ierr==9) then
703 
  call valdef_dc(sym,val)          ! Assign val to symbol sym
 703 
  call valdef_dc(sym,val)          ! Assign val to symbol sym
704 
end if
 704 
end if
705 
 
 705 
 
706 
end subroutine strdef
 706 
end subroutine strdef
707 
 
 707 
 
708 
 
 708 
 
709 
!**********************************************************************
 709 
!**********************************************************************
710 
 
 710 
 
711 
SUBROUTINE VALUEP(xinchar,cval)  ! Finds double precision complex value
711 
SUBROUTINE VALUEP(xinchar,cval)  ! Finds double precision complex value
712 
                                 ! corresponding to number string xinchar
712 
                                 ! corresponding to number string xinchar
713 
                                 ! or value in symbol table corresponding
713 
                                 ! or value in symbol table corresponding
714 
                                 ! to symbol name xinchar.
 714 
                                 ! to symbol name xinchar.
715 
 
 715 
 
716 
character (len=*):: xinchar
 716 
character (len=*):: xinchar
717 
complex(kc8) :: cval
 717 
complex(kc8) :: cval
718 
real(kr8) :: rval
 718 
real(kr8) :: rval
719 
 
 719 
 
720 
ierr=0
 720 
ierr=0
721 
 
 721 
 
722 
if(is_letter(xinchar(1:1))==.true.) then   ! xinchar is a symbol
 722 
if(is_letter(xinchar(1:1)).eqv..true.) then   ! xinchar is a symbol
723 
  call getparam(xinchar,cval)
 723 
  call getparam(xinchar,cval)
724 
else                               ! xinchar is a number string
 724 
else                               ! xinchar is a number string
725 
  call value(xinchar,rval,ios)     ! rval is the value of xinchar
 725 
  call value(xinchar,rval,ios)     ! rval is the value of xinchar
726 
  if(ios > 0) then
 726 
  if(ios > 0) then
727 
    ierr=3
 727 
    ierr=3
728 
    write(*,*) 'Error: number string ',trim(xinchar),' does not correspond to a valid number'
728 
    write(*,*) 'Error: number string ',trim(xinchar),' does not correspond to a valid number'
729 
    write(*,*)
 729 
    write(*,*)
730 
  end if
 730 
  end if
731 
  cval=cmplx(rval,0.0_kr8,kc8)
 731 
  cval=cmplx(rval,0.0_kr8,kc8)
732 
  return
 732 
  return
733 
end if
 733 
end if
734 
 
 734 
 
735 
end subroutine valuep
 735 
end subroutine valuep
736 
 
 736 
 
737 
 
 737 
 
738 
!**********************************************************************
 738 
!**********************************************************************
739 
 
 739 
 
740 
 
 740 
 
741 
SUBROUTINE PUSHOP(op)  ! Puts an operator on operator stack
 741 
SUBROUTINE PUSHOP(op)  ! Puts an operator on operator stack
742 
 
 742 
 
743 
type(item):: op
 743 
type(item):: op
744 
 
 744 
 
745 
itop=itop+1
 745 
itop=itop+1
746 
if(itop > numtok) then
 746 
if(itop > numtok) then
747 
  write(*,*) 'Error: operator stack overflow in evaluation of expression'
 747 
  write(*,*) 'Error: operator stack overflow in evaluation of expression'
748 
  write(*,*)
 748 
  write(*,*)
749 
  return
 749 
  return
750 
end if
 750 
end if
751 
opstack(itop)=op
 751 
opstack(itop)=op
752 
 
 752 
 
753 
end subroutine pushop
 753 
end subroutine pushop
754 
 
 754 
 
755 
SUBROUTINE POPOP(op) ! Takes top operator of operator stack and assigns it to op
 755 
SUBROUTINE POPOP(op) ! Takes top operator of operator stack and assigns it to op
756 
 
 756 
 
757 
type(item):: op
 757 
type(item):: op
758 
 
 758 
 
759 
op=opstack(itop)
 759 
op=opstack(itop)
760 
itop=itop-1
 760 
itop=itop-1
761 
 
 761 
 
762 
end subroutine popop
 762 
end subroutine popop
763 
 
 763 
 
764 
SUBROUTINE PUSHVAL(val) ! Puts value on value stack
 764 
SUBROUTINE PUSHVAL(val) ! Puts value on value stack
765 
 
 765 
 
766 
complex(kc8) :: val
 766 
complex(kc8) :: val
767 
 
 767 
 
768 
itop=itop+1
 768 
itop=itop+1
769 
if(itop > numtok) then
 769 
if(itop > numtok) then
770 
  write(*,*) 'Error: value stack overflow in evaluation of expression'
 770 
  write(*,*) 'Error: value stack overflow in evaluation of expression'
771 
  write(*,*)
 771 
  write(*,*)
772 
  return
 772 
  return
773 
end if
 773 
end if
774 
valstack(itop)=val
 774 
valstack(itop)=val
775 
 
 775 
 
776 
end subroutine pushval
 776 
end subroutine pushval
777 
 
 777 
 
778 
SUBROUTINE POPVAL(val) ! Takes top value off value stack and assigns it to val
 778 
SUBROUTINE POPVAL(val) ! Takes top value off value stack and assigns it to val
779 
 
 779 
 
780 
complex(kc8) :: val
 780 
complex(kc8) :: val
781 
 
 781 
 
782 
val=valstack(itop)
 782 
val=valstack(itop)
783 
itop=itop-1
 783 
itop=itop-1
784 
 
 784 
 
785 
end subroutine popval
 785 
end subroutine popval
786 
 
 786 
 
787 
!**********************************************************************
 787 
!**********************************************************************
788 
 
 788 
 
789 
SUBROUTINE GETPARAM_DC(sym,var)  ! Find double precision complex value var
 789 
SUBROUTINE GETPARAM_DC(sym,var)  ! Find double precision complex value var
790 
                                 ! corresponding to symbol sym
 790 
                                 ! corresponding to symbol sym
791 
 
 791 
 
792 
character(len=*) :: sym
 792 
character(len=*) :: sym
793 
character(len=len_trim(sym)) :: usym
 793 
character(len=len_trim(sym)) :: usym
794 
complex(kc8) :: var
 794 
complex(kc8) :: var
795 
 
 795 
 
796 
ierr=0
 796 
ierr=0
797 
sym=adjustl(sym)
 797 
sym=adjustl(sym)
798 
if(is_letter(sym(1:1))==.false. .or. len_trim(sym)>24) then
 798 
if(is_letter(sym(1:1)).eqv..false. .or. len_trim(sym).gt.24) then
799 
  ierr=11
 799 
  ierr=11
800 
  write(*,*) 'Error: symbol ',trim(sym),' has incorrect format'
 800 
  write(*,*) 'Error: symbol ',trim(sym),' has incorrect format'
801 
  write(*,*)
 801 
  write(*,*)
802 
  return
 802 
  return
803 
end if
 803 
end if
804 
ifind=0
 804 
ifind=0
805 
usym=uppercase(sym)
 805 
usym=uppercase(sym)
806 
do j=1,nparams
 806 
do j=1,nparams
807 
  if(trim(usym) == trim(params(j)%symbol)) then
 807 
  if(trim(usym) == trim(params(j)%symbol)) then
808 
    var=params(j)%value
 808 
    var=params(j)%value
809 
    ifind=j
 809 
    ifind=j
810 
    exit
 810 
    exit
811 
  end if
 811 
  end if
812 
end do
 812 
end do
813 
if(ifind == 0) then
813 
if(ifind == 0) then
814 
  ierr=4
 814 
  ierr=4
815 
  write(*,*) 'Error: symbol ',trim(sym), ' not in symbol table'
 815 
  write(*,*) 'Error: symbol ',trim(sym), ' not in symbol table'
816 
  write(*,*)
816 
  write(*,*)
817 
  return
 817 
  return
818 
end if
 818 
end if
819 
 
 819 
 
820 
end subroutine getparam_dc
 820 
end subroutine getparam_dc
821 
 
 821 
 
822 
!**********************************************************************
 822 
!**********************************************************************
823 
 
 823 
 
824 
SUBROUTINE GETPARAM_SC(sym,var)  ! Find single precision complex value var
 824 
SUBROUTINE GETPARAM_SC(sym,var)  ! Find single precision complex value var
825 
                                 ! corresponding to symbol sym
 825 
                                 ! corresponding to symbol sym
826 
 
 826 
 
827 
 
 827 
 
828 
character(len=*) :: sym
 828 
character(len=*) :: sym
829 
complex(kc4) :: var
 829 
complex(kc4) :: var
830 
complex(kc8) :: vard
 830 
complex(kc8) :: vard
831 
 
 831 
 
832 
call getparam_dc(sym,vard)
 832 
call getparam_dc(sym,vard)
833 
var=vard
 833 
var=vard
834 
 
 834 
 
835 
end subroutine getparam_sc
 835 
end subroutine getparam_sc
836 
 
 836 
 
837 
!**********************************************************************
 837 
!**********************************************************************
838 
 
 838 
 
839 
SUBROUTINE GETPARAM_DR(sym,var)  ! Find double precision real value var
 839 
SUBROUTINE GETPARAM_DR(sym,var)  ! Find double precision real value var
840 
                                 ! corresponding to symbol sym
 840 
                                 ! corresponding to symbol sym
841 
 
 841 
 
842 
 
 842 
 
843 
character(len=*) :: sym
 843 
character(len=*) :: sym
844 
real(kr8) :: var
 844 
real(kr8) :: var
845 
complex(kc8) :: vard
 845 
complex(kc8) :: vard
846 
 
 846 
 
847 
call getparam_dc(sym,vard)
 847 
call getparam_dc(sym,vard)
848 
var=real(vard,kr8)
 848 
var=real(vard,kr8)
849 
 
 849 
 
850 
end subroutine getparam_dr
 850 
end subroutine getparam_dr
851 
 
 851 
 
852 
!**********************************************************************
 852 
!**********************************************************************
853 
 
 853 
 
854 
SUBROUTINE GETPARAM_SR(sym,var)  ! Find single precision real value var
 854 
SUBROUTINE GETPARAM_SR(sym,var)  ! Find single precision real value var
855 
                                 ! corresponding to symbol sym
 855 
                                 ! corresponding to symbol sym
856 
 
 856 
 
857 
 
 857 
 
858 
character(len=*) :: sym
 858 
character(len=*) :: sym
859 
real(kr4) :: var
 859 
real(kr4) :: var
860 
complex(kc8) :: vard
 860 
complex(kc8) :: vard
861 
 
 861 
 
862 
call getparam_dc(sym,vard)
 862 
call getparam_dc(sym,vard)
863 
var=real(vard)
 863 
var=real(vard)
864 
 
 864 
 
865 
end subroutine getparam_sr
 865 
end subroutine getparam_sr
866 
 
 866 
 
867 
!**********************************************************************
 867 
!**********************************************************************
868 
 
 868 
 
869 
SUBROUTINE GETPARAM_DI(sym,ivar)  ! Find double precision integer value ivar
 869 
SUBROUTINE GETPARAM_DI(sym,ivar)  ! Find double precision integer value ivar
870 
                                  ! corresponding to symbol sym
 870 
                                  ! corresponding to symbol sym
871 
 
 871 
 
872 
 
 872 
 
873 
character(len=*) :: sym
 873 
character(len=*) :: sym
874 
integer(ki8) :: ivar
 874 
integer(ki8) :: ivar
875 
complex(kc8) :: vard
 875 
complex(kc8) :: vard
876 
 
 876 
 
877 
call getparam_dc(sym,vard)
 877 
call getparam_dc(sym,vard)
878 
ivar=nint(real(vard,kr8),ki8)
 878 
ivar=nint(real(vard,kr8),ki8)
879 
 
 879 
 
880 
end subroutine getparam_di
 880 
end subroutine getparam_di
881 
 
 881 
 
882 
!**********************************************************************
 882 
!**********************************************************************
883 
 
 883 
 
884 
SUBROUTINE GETPARAM_SI(sym,ivar)  ! Find single precision integer value ivar
 884 
SUBROUTINE GETPARAM_SI(sym,ivar)  ! Find single precision integer value ivar
885 
                                  ! corresponding to symbol sym
 885 
                                  ! corresponding to symbol sym
886 
 
 886 
 
887 
 
 887 
 
888 
character(len=*) :: sym
 888 
character(len=*) :: sym
889 
integer(ki4) :: ivar
 889 
integer(ki4) :: ivar
890 
complex(kc8) :: vard
 890 
complex(kc8) :: vard
891 
 
 891 
 
892 
call getparam_dc(sym,vard)
 892 
call getparam_dc(sym,vard)
893 
ivar=nint(real(vard,kr8),ki4)
 893 
ivar=nint(real(vard,kr8),ki4)
894 
 
 894 
 
895 
end subroutine getparam_si
 895 
end subroutine getparam_si
896 
 
 896 
 
897 
!**********************************************************************
 897 
!**********************************************************************
898 
 
 898 
 
899 
SUBROUTINE EVALEQN(eqn)  ! Evaluate an equation
 899 
SUBROUTINE EVALEQN(eqn)  ! Evaluate an equation
900 
 
 900 
 
901 
character(len=*) :: eqn
 901 
character(len=*) :: eqn
902 
character(len=len(eqn)) :: args(2)
 902 
character(len=len(eqn)) :: args(2)
903 
complex(kc8) :: val
 903 
complex(kc8) :: val
904 
 
 904 
 
905 
call parse(eqn,'=',args,nargs)   ! Seperate right- and left-hand-sides
 905 
call parse(eqn,'=',args,nargs)   ! Seperate right- and left-hand-sides
906 
call defparam(adjustl(args(1)),args(2)) ! Evaluate right-hand-side and
 906 
call defparam(adjustl(args(1)),args(2)) ! Evaluate right-hand-side and
907 
                                        ! assign to symbol on the
 907 
                                        ! assign to symbol on the
908 
                                        ! left-hand-side.
 908 
                                        ! left-hand-side.
909 
end subroutine evaleqn
 909 
end subroutine evaleqn
910 
 
 910 
 
911 
!**********************************************************************
 911 
!**********************************************************************
912 
 
 912 
 
913 
SUBROUTINE LISTVAR      ! List all variables and their values
 913 
SUBROUTINE LISTVAR      ! List all variables and their values
914 
 
 914 
 
915 
write(*,'(/a)') ' VARIABLE LIST:'
 915 
write(*,'(/a)') ' VARIABLE LIST:'
916 
if(nparams == 0) then            ! Initialize symbol table
 916 
if(nparams == 0) then            ! Initialize symbol table
917 
  params(1)%symbol='PI'
 917 
  params(1)%symbol='PI'
918 
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
 918 
  params(1)%value=(3.14159265358979_kr8,0.0_kr8)
919 
  params(2)%symbol='I'
 919 
  params(2)%symbol='I'
920 
  params(2)%value=(0.0_kr8,1.0_kr8)
 920 
  params(2)%value=(0.0_kr8,1.0_kr8)
921 
  nparams=2
 921 
  nparams=2
922 
end if
 922 
end if
923 
do i=1,nparams
 923 
do i=1,nparams
924 
  write(*,*) trim(params(i)%symbol),' = ',params(i)%value
 924 
  write(*,*) trim(params(i)%symbol),' = ',params(i)%value
925 
end do
 925 
end do
926 
 
 926 
 
927 
end subroutine listvar
 927 
end subroutine listvar
928 
 
 928 
 
929 
!**********************************************************************
 929 
!**********************************************************************
930 
 
 930 
 
931 
end module evaluate
 931 
end module evaluate