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