MSU

LIBRARIES
“

 

 

RETURNING MATERIALS:
Place in book drop to
remove this checkout from
your record. FINES wiII
be charged if book is
returned after the date
stamped beIow.

 

 

 

 

 

IMPLEMENTING A FUNCTION LIBRARY
FOR
NONLINEAR DYNAMIC SYSTEMS

BY

Tsun-Yu Chou

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Mechanical Engineering

1983

IJI'V"

ABSTRACT

IMPLEMENTING A FUNCTION LIBRARY
FOR
NONLINEAR DYNAMIC SYSTEMS

BY

Tsun-Yu Chou

ENPORT-S is a bond-graph based interactive simulation
program with extensive graphic ability that allows the user
to control and analyze the simulation of physical systems
effectively. To extend this ability to the nonlinear domain
an immediate need is an explicit expression of the nonlinear
constitutive equations that can be defined by the user and
modified easily and efficiently. An interactive program,
NOLMOD, has been written which is based on the core concept
of a prestored function library. Certain efficiencies are
available by using a pre-programmed library; flexibility

must be achieved by other special methods.

Chapter

TABLE OF CONTENTS

LIST OF TABLES o o o o o o o o o o o o 'o o o o o o Io

LI ST OF FIGURES O O O O O O O O O O O O O O O O O O O

I NTRODUCT I ON I O O O O O O I O O O O O O O O O

EQUATION FORMULATION FOR NONLINEAR BOND GRAPHS

2.1
2.2
2.3
2.4
2.5

Bond graph structure . . . . . . . . . .
Key variables and causal considerations
System equations and state equations . .
Example of matrix formulation . . . . .

Computing implications . . . . . . . . .

EXAMPLES O O O O O O O I O O O O O O O O O O O

3.1
3.2
3.3

Spring-mass-damper model . . . . . . . .
Hydraulic system . . . . . . . . . . .g.

Dry friction O O O O O O O O O O O O O O

DES I GN OF PROGRAM NOLMOD O C O O O O C I O O O

4.1
4.2
4.3
4.4
4.5

Basic considerations . . . . . . . . . .
Classification of equations . . . . . .
Function library . . . . . . . . . . . .
Evaluation of functions. . . . . . . . .

Structure of NOLMOD . . . . . . . . . .

ii

Page

iv

(DONUIU'IH

12
17
18
18
23
25
28
28
29
31
33
34

Chapter . Page

5. CONCLUSION . . . . . . . . . . . . . . . . . . . 39

5.1 Conclusion . . . . . . . . . . . . . . . . 39

5.2 Limitations and future work . . . . . . . 40

LIST OF REFERENCES . . . . . . . . . . . . . . . . . . 42
APPENDICES

A. HOW TO USE NOLMOD . . . . . . . . . . . . . . . Al

B. IMPLEMENTATION OF A NEW FUNCTION . . . . . . . Bl

C O S OURC E CODE FOR NOLMOD O O O O O O O O O O O O O C 1

iii

LIST OF TABLES

Table Page
1. List of library functions . . . . . . . . . . . . 32

2. List of functions in alphabetic order . . . . . . A13

iv

Figure

2.1
2.2

2.3

2.4

LIST OF FIGURES

Basic fields of a multiport system . . .

Identification of key vectors in a mixed
causality system . . . . . . . . . . . .

A multiport system having only integral
causality O O O O I O O O O O O O O O O

(a) Spring-mass-damper model
(b) Bond graph model . . . . . . . . . .

Augmented bond graph . . . . . . . . . .
Subroutine FCT . . . . . . . . . . . . .
Displacement diagram . . . . . . . . . .
Velocity diagram . . . . . . . . . . . .
Hydraulic system . . . . . . . . . . . .
Bond graph of hydraulic system . . . . .

(a) Bond graph of bent hose
(b) Modified bond graph of hose . . . .

Dynamic model of a vehicle . . . . . . .

(a) Bond graph of vehicle model
(b) Saturation function . . . . . . . .

(a) Modified bond graph of vehicle
(b) Inverse of saturation. . . . . . . .

Page

11

13
13
20
21
22
24
24

24
26

26

26

Figure Page
4.1 Basic components of NOLMOD . . . . . . . . . . . 35
4.2 Subroutines calling structure for NOLMOD . . . . 36
4.3 Data flow for function library . . . . . . . . . 38
A.l Function A352 . . . . . . . . . . . . . . . . . All
A.2 Function ABSQRT . . . . . . . . . . . . . . . . All
A.3 Function RCPSQR . . . . . . . . . . . . . . . . All
A.4 Function COULOM . . . . . . . . . . . . . . . . A12
A.5 Function BRKPT . . . . . . . . . . . . . . . . . A12
A.6 Function SATUR . . . . . . . . . . . . . . . . . A12
A.7 Function ADJSAT . . . . . . . . . . . . . . . . A12
A.8 Function PWLIN . . . . . . . . . . . . . . . . . A12

A09 Function INTPLO O O O O O C O O O O O O O O O 0 A12

vi

CHAPTER 1

INTRODUCTION

During the past decade the digital computer has come to
play a very important role in system simulation. Many
general purpose as well as special purpose modeling and
simulation programs have been developed. One general
requirement of such a program is that the user need not have
extensive knowledge of computer programming and numerical
techniques. The typical user can learn to operate such a
program, if it is well-designed, within a very short period

of time.

Many important engineering and physical systems are
described by differential equations. Simulation programs
provide a way of solving such differential equations and
also provide graphical and tabular output that makes it very
easy for the user to study the behavior of the systems.
Selected parameters and functions can be entered to provide
necessary information about the differential equations.
They are constructed to define the relationships between

dependent and independent variables and among dependent

variables. How to provide such information to a computer is
a topic in which we are interested. Therefore several
simulation programs are studied, including CSMP,
SUPER-SCEPTRE and DIFFEQ. CSMP is a general purpose
block-diagram-oriented program developed by IBM for
simulation of physical and engineering systems '[1-3].
SUPER-SCEPTRE is a special purpose program intended
primarily for electrical systems[4]. DIFFEQ is a general
equation-based program developed at Michigan State

University for solving ordinary differential equations [5].

Patterns for entering data are very similar in both
CSMP and SUPER-SCEPTRE. Control statements and data
statements are used as input data. Control statements tell
the system what kind of data is provided and how to handle
it. Data statements give the value of parameters, the types
of functions and other related information. Some commonly
used functions are pre-programmed, and both programs also
provide function-generating routines for special functions.
Then the system, based on the data given, prepares the
necessary subroutines for calculations. The user can also
write his own subroutine, if he is familiar with computer
programming. Then the program will compile and load the
binary file and set up the run-time program. This is the
approach taken in DIFFEQ. A FORTRAN subroutine is written
describing the differential equations and the system

provides the calculation and display capabilities.

The common point of all three programs is that they
have to compile the function subroutines and combine them
with the rest of the program. There are several
disadvanteges of such a process. First of all, it take both
time and disk space to do it. When a lot of users are
active at one time and every user stores a binary file of
the program on the disk, that can use up a lot of space.
And the compiling and loading process takes computer time.
An advantage, however, is that the run-time program can be
saved for re-run with new data and execution control
statements. Subsequent simulation jobs using the same model
may avoid repeating the translating, compiling and linking
process. Some careful planning is necessary, coupled with a
sound background in FORTRAN, to ensure that the saved model
has sufficient capacity for later uses. Since both CSMP and
SUPER*SCEPTRE are used in batch mode, user has to submit the
job all over again when the data are changed. But DIFFEQ is
a fully interactive program in which changing of parameters
can be done on-line. This feature enables the user to have
direct interaction with the simulation. This is a great
advantage in design, since it can be re-run relatively
quickly until the user is satisfied with the results. But
changing of a function type in DIFFEQ must be done by
modifying the subroutine. Then the compiling process must

be repeated.

ENPORT is an interactive program based on bond graph

theory for modeling and simulation of physical systems which
has many of the desirable properties mentioned above [6-7].
However, at this time the ENPORT program can only handle
linear or linearized models. Demand that ENPORT be able to
handle nonlinear models is increasing. The purpose of this
work is to develop a program that extends ENPORT into the
nonlinear domain. For nonlinear problems functions are
specified, instead of just constant parameters. Therefore
our strategy is to set up a function library for such a
purpose. The program NOLMOD was written to achieve this
goal. The user can choose desired functions from the
library and can modify these choices without leaving the
program. To use this system the designer first runs ENPORT
to read the bond graph as original data and to assign the
causality. Then NOLMOD takes over the rest of the work.
The program NOLMOD has been installed on the PRIME-750
computer of A. H. Case Center, College of Engineering,

Michigan State University.

CHAPTER 2
EQUATION FORMULATION FOR NONLINEAR BOND GRAPHS
2.1 Bond Graph structure

One of the most important features of the. bond graph
method is the concept of fields and junction structures.
These terms were defined by H. Paynter [8]. A multiport
system is partitioned into a set of distinct, interconnected
fields, namely, a source field, a storage field, a
dissipation field, and a junction structure. By definition
a junction structure is a collection of junction elements :
0- and 1- junctions, and transformers (TF) and gyrators
(GY). In a standard bond graph every bond is connected to a
junction element at least at one end. It is useful to
classify the graph bonds into external bonds and internal
bonds. A bond connected between a field element (C, I, R,
SE, SF) and a junction element (0, 1, TF, GY) is called an
external bond. An internal bond is a bond connected between
two elements of the junction set. Furthermore, external
bonds are classified into three groups according to the
nature of the elements they adjoin : source field bonds (SE

and SF), storage field bonds (C and I) and dissipation field

bonds (R). The equation structure implied by a bond graph
is constructed based on these classifications. Figure 2.1

shows the bond graph field structure [6].

2.2 Key variables and causal considerations

Figure 2.2 shows the key variables (vectors) used for
representing the input and output variables for each field.
U is the input to the junction structure from the source
field, V is the output from the junction structure to the
source field; these are associated with the source field
bonds. 2 is the co-energy variable vector, X dot is the
time derivative of the energy variables vector; these are
associated with the storage field bonds. Do is the output
from the dissipation field to the junction structure, Di is
the input to the dissipation field from the junction
structure; these are associated with the dissipation field
bonds. Subscripts i and d stand for integral and derivative

ports of the storage field, respectively.

Causal consideration is a very important aspect in the
bond graph method. It provides all the information
concerning input and output variables. Some multiports are
strongly constrained to certain causality (SE and SF). For
C- and I-fields, the constitutive laws are quite different
for different causalities. Usually integral causalities are

sought. Derivative causality arises for systems in which

Storage
field

 

C

 

Source field

I

 

 

 

 

 

 

 

 

 

 

 

 

(1, O, TF, GY‘ °

 

 

SE SF
...... Dissipation
field
Junction R
Structure

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 2.1 Basic fields of a multiport system.
Source field
(SE, SF)
1
U V
‘ i
i Junction d
(0.1) Structure (C,I)
. Z
1+(1. 0. TF, GY) ‘1
Integral D I D Derivative
storage i‘ 0 storage
Dissipation
field (R)

 

 

 

Figure 2.2 Identification of key vectors in a mixed
causality system.

some storage elements are not dynamically independent. The
energy variable of the derivative storage element port will
be expressed algebraically in terms of other independent
energy variables. Therefore any storage element ports with
derivative causality do not contribute to the state variable
set. For a nonlinear system it may be difficult to express
one energy variable in terms of the others explicitly, since
they might be coupled by nonlinear algebraic equations.
This situation is typical of nonlinear mechanics, for
example. The causality for the R-field usually is
determined by the source and storage-field elements.
Usually they are relatively indifferent to causality.
However, in the nonlinear case some constitutive laws for
the R-field are not uniquely determined for certain

causalities. An example will be shown in Chapter 3.

2.3 System equations and state equations

First we shall look at the situation for linear
systems, useful as a guide to our thinking in the nonlinear
case. The field equations for linear systems are as

follows:

Source field :
U s U(t) (2.1)
Storage field :

Z=FX+FX (2.2a)
i ii 1 id 6

and

Z = F X + F X (2.2b)
d di 1 dd d

Dissipation field :

o . L D (2.3)
o i

Junction structure equations :

x = s z + s x + s D + s u (2.4a)
i 11 i 12 d 13 o 14

z = s z + s i + s D + s U (2.4b)
d 21 i 22 d 23 o 24

D = s z + s x + s D + s U (2.4c)
1 31 1 32 d 33 o 34

v = s z. + s x + s D + s U (2.4a)

The system equations above can be reduced to a single
matrix state-space equation by some suitable manipulations.
(We assume the required inverses exist.) Thus we arrive at

x. = A x. + B U + E U (2.5)
l 1

The E matrix only appear when we have mixed causality.
One may expect that all entries in matrices A, B and E are
_real numbers which are derived from the parameters of the
physical elements. Matrices A, B and E can be calculated

from a causal bond graph after parameters are given. With

10

given X and U, X dot can be obtained easily. And by use of
the equations (2.1), (2.2), ~(2.3) and (2.4), all system

variables can be found.

For some nonlinear systems equations (2.2) and (2.3)

are no longer valid. They will be in following form:

For the storage fields:

2 = (X ,X ) (2.6a)
i ¢Fi i d

z = 4’ (x ,x) (2.6b)
d Fd i d

For the dissipation field:

D =<F (D ) (2.7)
o L i
In this case state-space equations of the form:

in. =<i>(x.,u) (2.8)
1 1

may not be obtainable. System equations may be firmly
coupled together. It may be very complex to reduce the'
system equations to a single state-space equation in the

form of equation (2.8). Typical stumbling blocks are

equations (2.6) and (2.7).

Therefore in this work we restrict the range of

consideration to the case of:

(l) integral causality (hence xg=0)

(2)

(3)

See

11

no implicit R-fields (hence S =0)
33

constant moduli for TF and CY (hence S have constant
13
elements)

Figure 2.3 for the key vectors in this case. Then the

 

Source field

 

 

 

 

 

 

 

 

 

 

 

 

(SE, SF )
I
Storage U V Dissipation
field , , field
(C I X Junction Di (R)
' Structure
Z DO
‘ (1, O, TF, GY ~

 

 

 

 

 

 

 

 

 

Figure 2.3 A multiport system having only integral

causality.

field equations and junction structure equations will be

simplified as follows:

2 = 4: (x) (2.9)
F
D a ct (D ) (2.7)
o L i
i = s z + s D + s U (2.10a)
11 13 o 14
D = s z + s U (2.10b)

12

Rearrange the equations in the following sequence:

U = U(t) (2.1)
z = 4>(x) (2.9)
F
D = s z + s U (2.10b)
i 31 34
D = (D ) (2.7)
o #L i
i = s z + s D + s U (2.10a)
11 13 o 14

With U and X given, then 2, Di, Do and X dot can be
calculated step by step according to the above sequence
[11]. Viewed as a total procedure this is equivalent to the

explicit form given by equation (2.8).

2.4 Example of matrix formulation

An example is shown in this section to demonstrate the
reduction procedures for both linear and nonlinear models.
Figure 2.4 is a simple example of a spring-mass-damper
model. An augmented bond graph is shown in Figure 2.5 which
displays the fields and junction structure. The field

equations for linear system are as follows:

Source field :

U = U(t) (2.1)

 

 

 

 

 

 

 

 

Figure 2.4a

13

Figure 2.4b

Figure 2.4 (a) Spring-mass-damper model
(b) Bond graph model.

 

   

 
 

 

 

I u
i Se : Source field
i ___________ J

F"'1

l I 2 1

: C

I I f ----------- 1

I I I

i i i J! { Junction

E E E 1 E Structure

5 M)

i I u

L___J 3 4

Storage { ----- 3----1

field I R u Dissipation
2 : field
h- --------- 4

Figure 2.5 Augmented bond graph model.

14

Storage field :

z = F x

F k 0 q
2 __ 2

v o l/m p
3 3

Dissipation field :

[2% PM]

 

 

 

 

Junction structure equations :
3 Is 5
== ll 13
D S S
i _ 31 33
q I0 1 g o g
2 I l
| I
P = '1 0 :"1 :
3 ........ J----J
-.--_ . .
V 0 l : 0 :
. 41
The system equations can be

equation:

 

reduced

 

 

 

 

to a

(2.11a)

(2.11b)

(2.12a)

(2.12b)

(2.13a)

(2.13b)

state-space

15

x. = A x + B U (2.14)

Where :

-1
A = (:5 + s L ( I - s L) s ‘]F (2.15a)
11 13 33 31

-1
B = [s + s L ( I - s L) s ] F (2.15b)
14 13 33 34

Therefore, by matrix operations on equations 2.11a, 2.12b

and l.3b,
[ 0 l/m
A = (5.16a)
-k -b/m
0
B = [ J (5.16b)
1

If the storage and dissipation elements are nonlinear,
then an explicit form of state-space equation (2.14) may not
be obtainable. We must follow the procedures described in
section 2.3 to obtain the implicit state-space equation.

The field equations will be as follows:
Source field :

U = U(t) (2.1)
Storage field :

z = d}(x) ‘ (2.9)

16

F = (q ) (2.17a)
2 4E 2

V = (P ) (2.17b)
3 43 3

Junction structure equations :

D = S 2 + S U (2.10b)
i 31 34 '

D = [o 1 ] .F ‘ + [o ] U (2.18)

 

 

Dissipation field :

D = 4>(D ) (2.7)
o L i

F = (v ) (2.19)
44:.

Implicit state-space equation :

 

 

 

 

 

 

 

 

x = s z + s D + s U (2.10a)
11 13 o 14
[.4 I 0 1‘ PF 7 r o. [F ] i 0 [U]
2 2 4
= + + (2.20)
e -1 o v -1j _ 1}
3] ' 1 ~ 3‘ '

 

 

Values of F2, V3 and F4 should be obtained from equations

(2.17), (2.18) and (2.19).

 

_

 

17

2.5 Computing implications

We now turn our attention to the practical question of
how to implement equations (2.7) and (2.9) and the procedure
for solving them. For equation (2.7) and (2.9), the
nonlinear functions describe the relations between input and
output variables on a port-by-port basis. In order for the
engineer to supply such information to the program
interactively, functions may be pre-programmed and
parameter-type data given during interactive running. The
most common engineering models involve one-port elements, in
which case only one dependent variable is involved. For
multi-port field elements the output set is a function of
the input set of two or more variables. It is more
difficult to define a set of library functions that can be
used for describing a variety of multi-port elements.
Therefore we restrict the bond graph to contain one-port

field elements only. This will simplify the problem again.

A function library that contains commonly used
functions and a function-generating routine are used to
serve the purpose. The library can also be used to describe
nonlinear transformer and gyrator element moduli, although

they are not implemented at this stage of the work.

CHAPTER 3

EXAMPLES

3.1 Spring-mass-damper model

For better understanding of the procedures for
simulating nonlinear systems by the bond graph method and
NOLMOD program, the example of a spring-mass-damper system
in Section 2.5 is used for demonstration. Commands and data

'listed below are given to ENPORT.

GRAPH

SE 1 , C 2 , I 3 , R 4 , 1 1 2 3 4 .
CAUSAL

After causality has been assigned, NOLMOD will classify the
basic equations according to types of elements and causal

conditions.
INTEGRAL STORAGE ELEMENTS
82=FCN2(Q2)
F3=FCN3(P3)

DISSIPATION ELEMENTS

18

19
E 4 = FCN 4 (F 4)

SOURCE ELEMENTS
E 1 = FCN ( TIME )
Then the user can define the functions for the storage and
dissipation constitutive equations. The user-defined

functions are as follow :

FUNCTION 2 :

E 2 = 0.0000E-01 + 4.0000E+00 * Q 2

FUNCTION 3 :

P 3 = 0.0000E-01 + 6.0000E-01 * F 3

FUNCTION 4 :

E 4 = 8.0000E-01 * s:ou( F 4 ) * SQRT( ABS( F 4 ) )

Following the causal conditions in Figure 4.5 the state

equations will be :
o = F (3.1)

P = E r E - E (3.2)
3 1 2 4

Then we can integrate equations (3.1) and (3.2) to obtain
the displacement and momentum of the system. Now we may use
DIFFEQ to simulate the system for demonstration. Figure 3.l
is the subroutine FCT used in DIFFEQ [4] and Figure 3.2 and

3.3 are time responses of the displacement and momentum.

20

.90h mcfipzomnsm H.m mhswwm

atu
tcahu- 0
cm I mu I «u o .m.>¢ua
nu o .«.>¢ua u
mzouhcaou wocnmIUhCPmIIIw
.nu.mo¢ a nu u a - vu u
. onrcnow adunu zanhcmnmmnaIuow
Amv> a At\.«. a nu
an.) n v a mu
0
. mzouk¢aow agunu uo¢¢ormnuuw
m.a u c a «u o
. onhcaou oamnu uocaomIIIw
Anvcaon
Amvcuu:
Anvca-x o

baazn «hcounuo
uuuuuuunuuununnuuuuunuuunuunuuuuunuuao

kzmnonuuuoo zonku~au - Ianu
mmc: - : IIIo

hrcrmzoo uzuaam - x IIIW

mumou zo~ku~au . 'u IIIo
sztuoaummua mmcc . nu IIIo

women oznunm . mm .100

mmct uo hxonw: - aw IIIW

0.0 camv> u
0.“ can.) . mzouhnozoo gauhuznuncw

rakzmto: mac: n nmv> IIIu
zomhuuuumo etuunm a an.) 14-0
-uuuunuuuuunuuunuuuuuuuuunnuu-uuuuuuo
caxsourou\10ttoo
t.x Acwc
Aomvcn new“) Remv>¢uu vuatuu

.>¢ua.>.x.kou uzuraocnam

21

am

mm

.Emnmmflu EDPCmEoz m.m mhzmwm

mz_k
sm m_ a“

_

 

 

CJ*“UDCL_J<(LJUU:ELUZZF—

22

am

mm

.Emmwmflv HCmEmomHQmHQ m.m omsmwm

wzah
sm ma 5H

 

EECDIEUUZZF—IDIE

23

3.2 Hydraulic system

Figure 3.1 shows a simple hydraulic system. A
hydraulic motor is driven by a pump. The 2-block 4-way
valve is in on position, a pressure relief valve is used to
maintain system pressure and an adjustable throttle is used.
to control the flow through the hydraulic motor. The
corresponding bond graph for such a hydraulic system is
shown in Figure 3.2. Causality is uniquely assigned
according to basic rules. In modeling the pipe line,
attention is focused on the bent hose, since the pressure
drop is significant at the bent compared to a straight line.
Figure 3.3a shows a portion of bond graph for the bent hose
which includes inertia, compliance and dissipation.

Pressure drop in a bent tube is expressed by :
Q 2
P = K J?- (---) (3.3)

Figure 3.3a shows the causality of the dissipation element
Rh which indicates flow rate Q is input and pressure P is
output. Equation 3.1 matches the causal condition of the R
element, namely, pressure P (left-hand side) is given

explicitly by flow rate Q (right-hand side).

If the inertia of the fluid is considered to be

insignificant in the system, one might like to exclude the

24

 

Figure 3.4 Hydraulic system.

Rp Rpf Rh Ch Rm CS R1
1 I I I, I I I I I
st_.1I—-TFI—~o-=IlI->o—-I0—->I1-=ITF—-ll)—>o—~I 1I—>Se
I L l l
Ip Ih Im Il
pump hose relief motor shaft load

valve

Figure 3.5 Bond graph of hydraulic system.

Rh Ch
I. l Rh Ch
'—=fl 1)-" 0'-=d 'r ‘I
i —»I 1I—-=-0——>(
I
h
Figure 3.6a Figure 3.6b

Figure 3.6 (a) Bond graph of bent hose
(b) Modified bond graph of hose.

25

inertia from the model. Figure 3.3b shows the modified bond
graph after the I element is removed. The causal condition
of the Rh element is changed, while the rest of the system
remains the same. Therefore the input and output relation
is changed. Equation 3.1 can be inverted easily to fulfill
the causal condition. For the user's convenience it is a
good practice to have the simulation program invert the
function when requested. The user can ask the program to
invert the input-output relation of the function, if
required, while supplying the same function as before.
However, for some cases the inverse of a function might not
exist. Therefore an alternative can be very helpful. An

example of this situation is presented in the next section.
3.3 Dry friction

A vehicle can be described by a simple model consisting
of masses, springs and a damper as shown in Figure 3.4.
Here M is the mass of the vehicle and m is the mass of the
tire. The spring function k models the tire. A bond graph
model of such a system is expressed in Figure 3.5a, in which
causalities are also assigned. Now if we assume the damper
is not well lubricated, such that it experiences dry
friction when the car is running. The relation between
force and velocity for dry friction can be approximately
expressed by the saturation function (Figure 3.5b). When

the velocity reaches a certain value the friction force will

 

 

 

 

 

Figure 3.7 Dynamic model of a vehicle.

Se-——(1-——(I
I /R
OF-wv1A\\-
l C F)

Se-filfill F, ......

T 4v

01—..0

I

S

 

 

f
Figure 3.8a Figure 3.8b

Figure 3.8 (a) Bond graph of vehicle model
(b) Saturation function.

 

 

Se
J. VI)
1-—=HI
CHI-fil/R A IF
I \c
Sf
Figure 3.9a Figure 3.96

Figure 3.9 (a) Modified bond graph of vehicle
(b) Inverse of saturation.

27

remain constant (Fs). According to the causality of the R
element, velocity is the input to the dissipation element
and force is output. The saturation function is suitable

for defining the relation between force and velocity.

Often the mass of the tire can be neglected compared to
the large mass of the vehicle. A modified bond graph model
is shown in Figure 3.6a. The causal condition for the R
element has changed. However, the inverse of the saturation
function has an undefined region for F, which is now the
input variable (Figure 3.6b). From the standpoint of
computation serious difficulties might arise in such a
system. Therefore an alternative might be necessary. One
way of doing that is to define a very large finite value for

the slope (instead of infinity).

The examples above show that causal conditions can
change when the physical model of the system is modified.
Nevertheless, we desire that the simulation program has the
ability to handle it. It should not only provide a variety
of functions but also be flexible in their causal

implementation.

CHAPTER 4

DESIGN OF PROGRAM NOLMOD

4.1 Basic considerations

NOLMOD, which can be imbedded in the ENPORT-5 program,
is written in standard FORTRAN V, which helps to make it
readily portable. The program has been designed to allow
easy implementation on a variety of computers. It has been
developed as a completely separate module, which can be used
as part of ENPORT-5 or used independently. When used with
ENPORT, a new command has been added for invoking NOLMOD
(see Appendix A for usage of program). Data are transmitted
through common blocks, which can be passed through a memory
or disk file interface between ENPORT and NOLMOD. Also it

is user-oriented for easy operation.

In ENPORT, bond graphs are taken as input data, then
power directions and causality are assigned. In linear
models, the next step is to enter parameters of constitutive
equations (i.e. for R, C and I, TF and CY). For the
nonlinear case, functional relationships between input and

output variables must be defined for the constitutive

28

29

relations. NOLMOD takes over the nonlinear part at this
point. Since we currently restrict the nonlinear elements
to be one-ports, only one input and one output variable is
required for each element. The basic relation is:
output = function (input)

The most important information for defining such a
relationship specifically is causality and the bond types
(or element types). Such information is supplied by ENPORT.
Based on such information, NOLMOD makes an explicit
classification of the basic equations. Both integral and
derivative causality, nonlinear 2-port -TF- and -GY-, and
multi-port elements are classified within the program. But
some classifications are reserved for future implementation.
Attention is focussed below on one-port integral storage

elements and dissipators.

4.2 Classification of equations

In the discussion below the following definitions are

used:
E = generalized effort,
F = generalized flow,
P = generalized momentum,
and Q = genaralized displacement.

According to the different types of elements and causal

conditions, equations are classified into the following

30

eight types:
1. Integral C elements :
E1 = FCN(Qi)
2. Integral I elements :
Pi = FCN(Pi)
3. Derivative C elements :
Qi = FCN(Ei)

4. Derivative I elements :

pi = FCN(Fi)
5. R elements :
Pi = FCN(Ei)
or
Ei = FCN(Fi)
6. Source elements :
Ei = FCN(Time)
or
Pi = FCN(Time)

7. Transformers :

81 = FCN(Ej), Fj = FCN(Fi)
or Pi = FCN(Fj), Ej = FCN(Ei)
8. Gyrators :
Pi = FCN(Ej), Fj = FCN(Ei)
or E1 = FCN(Fj), Ej = FCN(Fi)

where i and j are the particular bond indices.

In the current NOLMOD implementation source elements
(type 6) are handled separately, since usually they are
functions only of the indenpendent variable time. ' ENPORT

already has an efficient procedure for handling source

31

elements, therefore NOLMOD excludes them. For the other

types a function library has been built up.

Once the graph causality conditions are known, the
NOLMOD program classifies the system-formatted equations
into the types indicated in the last section. Input and
output variables are defined for every constitutive
equation. In some cases the user might have chosen a
different input/output orientation for a particular element.
Or the causal relation is changed due to modification of the
physical model, as we mentioned in Chapter 3. Therefore,
for the user's convenience, he (or she) can ask the program
to invert the input/output orientation. The program accepts
the user-defined equation form and inverts it to obtain the
necessary data for the system-formatted equations. For
functions that can not be inverted, the program will prompt
the user with an error message and ask for further

instructions.

4.3 Function library

Some commonly encountered functions are pre-programmed,
available for user selection in defining input/cutout
relations. A list of the library functions is shown in
Table 1. Details can be found in Appendix A. The user can
check the properties of the function to see if it meets the

required behavior of the model. When the user inverts the

LINEAR
POLY
SIN
ASIN
COS
ACOS
DSIN
DASIN
DCOS
DACOS
EXP
ALOG
ABSZ
ABSQRT
RCPSQR
COULOM

BRKPT

SATUR

32

Table 1. List of library functions

DEFINITION

*40<|< ’4'4 K: *4 *< *4 '4 *4 K2 *4 N *< *< K: *4 K: *4 K: '4

*-

\\\\

In

c1*x

c1*x + C2*x**2 + C3*X**3 + C4*X**4
SIN(C2*X+C3) + C4
ASIN(C2*X+C3) + C4
cos(c2*x+c3) + C4
ACOS(C2*X+C3) + C4
SIN(C2*X+C3) + C4
ASIN(C2*X+C3) + C4
COS(C2*x+C3) + C4
ACOS(C2*X+C3) + C4
EXP(C2*X+C3) + C4
ALOG(C2*X+C3) + C4

x * ABS(X)

SIGN(X) * SQRT(ABS(X))
c1 / (C2+x**2)

SIGN(X)
C1*(x-x1), x <=x1
C2*(x-x1), x1< x

, x <=x1
((Yl-Y2)/(Xl-X2))*(X-X1), x1<x<x2
, x2<=x

33

Table l. (cont'd)

TYPE DEFINITION
ADJSAT Y = Yl + C1*(X-Xl) , X <=x1
Y = Y1 + ((Yl-Y2)/(Xl-X2))*(X-Xl), x1<x<x2
Y = Y2 + C2*(X-X2) , XZ<=X
PWLIN PIECEWISE LINEAR
INTPLO INTERPOLATION POLYNOMIAL
input/output orientation, the user-defined function

parameters can be read into the program. NOLMOD will invert
the function and calculate the appropriate parameters for
system-formatted equation. All the parameters are stored in
the system-formatted form. Some function inverses are not
available. Therefore it is user‘s responsibility to select
other functions that are better suited to serve his (her)

purpose.
4.4 Evaluation of functions

A program is provided for evaluation of the defined
functions. Such evaluations are based on system-formatted
equations (equations 2.7 and 2.8 of Chapter 2), which are
determined by causality. Given the values of input
variables, the values of output variables are obtained.

This program will be used later to evaluate state-space

34

equations as described in Chapter 2. This is why the
evaluation is based on system-formatted equations, not

user-defined equations. (See Appendix A for its usage.)

4.5 Structure of NOLMOD

Figure 4.1 shows the basic components of the NOLMOD
program. The subroutines enclosed in solid blocks refer to
the major parts. Detailed calling sequences are shown in
Figure 4.2. All important data are transferred through
common blocks COMFILE and COMAREA, which are found in almost
every subroutine. COMFILE is the common block set for
ENPORT-5 and COMAREA is the common block set for NOLMOD. In
general only certain control data are transferred through
arguments of calling statements. To have an efficient way
for handling the data from the terminal, two subroutines
-CMREAD and CMNCPR- are used as the free-reader part of the
program. All the data are entered as string variables
through CMREAD. Then commands will go through CMNCPR for
mapping onto the command list. Constant data go through a

decoding process.

In order for NOLMOD to execute successfully certain
data must be supplied by ENPORT or by a saved file of
ENPORT. Those data include the number of elements (NEL),
number of bonds (NBD), element type list (IELLST), types of

bonds (IBT), pointer list for NBIMX (NPTR), and the list of

35

 

 

 

 

 

INIT
Initialization

 

 

 

DEFINE

Classify
equations

 

 

 
 

 

 

 

 

 

 

LOAD

Read bond
graph data
from saved
file

 

   

 

 

 

 

 

FCNDEF

Define and
modify
nonlinear
functions

 

 

 

 

SHOW

Display
functions

 

 

  
  
 

 

 

 

ENPORT
NLEVAL
Evaluation
____________ .
SWITZ I
Additional :
NOLMOD "91313131-":
_____________ '
COMMON INFO I
BLOCKS: :Display power!
.direction andI
COMFILE Icausality I
I _____________ I
COMAREA
{UTILITYj , ............ I
I HELPNEW :

'Display menu)
Iand library I
:functions J

Figure 4.1 Basic components of NOLMOD.

36

MAINC’--(*UTILITY*)'--CMREAD,CMNCPR,ERRMSG

L-INIT
--DEFINE
r-SWITZ

”-LOAD- 'INIT

“DEFINE
”-INFO
“-HELPNEW

--FCNDEF- ’FCNWTR

 

brNLEVAL-[-FCNLOOK

-CMNZER(ENPORT)

'EQRVS ---FCNWTR

 

 

’GETVAL---FMIDV

 

-FCNLIB-r-FCNSUBl-- --ENCON
--FCNSUBZ-~
. --GENPAR
. r-IORVSR
b-FCNSUBZZJ
b-PARTAK
--SUBWTR
--coonsx
L-1NTERp--—DTABLE
P-SHOW ----- FCNLOOK---FCNSUBl,FCNSUB2,...,FCNSUBZZ

Figure 4.2. Subroutine calling structure for NOLMOD

37

bonds incident to each element (NBIMX). (Details regarding
usage can be found in subroutines RFADBG and CLASS of the

ENPORT program.)

Before reading the bond graph data into NOLMOD, INIT
should be used to set the NOLMOD into neutral status. LOAD
reads the bond graph data from a saved file. It will
initialize both ENPORT and NOLMOD before reading the data.
Then subroutine DEFINE classifies the basic equations, as
described in section 4.2. Subroutine FCNDEF and subsequent
subroutines shown in Figure 4.2 perform operations of
setting up equations, providing inverse functions, modifying
the functions, and verifying the inputs. SHOW is used to
display functions in system-formatted form or in
user-defined form. And NLEVAL provides on-line evaluation
of the functions as defined by the user. A block diagram in
Figure 4.3 shows the structure of function library. The
interested reader can refer to the source code listed in

Appendix C for further details.

 

 

 

 

 

 

38

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

COMMAND
FCNDEF
A A
—-—-—*E Define (g, I)
functions T
I I
I l
I J
L---~ INVERSE
r"‘* CONTROL
l I
i
: I
< 111 l
I
SHOW ./ W I;
‘ Di Splay \ 7)
functions
NLEVAL 44 I;
"——“—5’ Evaluate K:; 7)
functions

 

 

 

 

 

.. ..... J ..... .

Lgeleeleiiee'

GETVAL

.J

 

COMMON BLOCKS

FUNCTION
LIBRARY

and

UTILITY

 

Figure 4.3 Data flow for function library.

 

CHAPTER 5

CONCLUSION

5.1 Conclusion

For the engineering user who wishes to simulate bond
graph models in an interactive environment the ENPORT
program has proven useful. However, the need to include
nonlinear models has led to the design of a library-based
approach for specifying nonlinear constitutive equations.
Implemented as both a separate program (NOLMOD) and as a
feature of ENPORT, it allows the user to choose from a
variety of function types and to set the parameters of the
functions in the causal form he (or she) prefers. NOLMOD
then inverts the function to match the system causal
requirements, if that is necessary. Cases where an inverse
does not exist are identified and the decision control is

returned to the user.

A manual for the use of NOLMOD has been developed. It
defines the function library, the command set, and gives
examples of use. Since the program is designed to be used
interactively, it is menu-based. On-line error diagnostics

39

40

are provided. Anyone familiar with bond graph modeling in
general, and ENPORT usage in particular, will find it easy

to use NOLMOD.

5.2 Limitations and future work

The functions collected in the library are not
sufficient to cover all the different needs that arise in
practice. It is expected that the function library will be
modified to meet additional needs from time to time. In the
future some special functions might be included to fulfill
the requirements of different engineering fields. Although
the program is well structured for ease of implementation, a
programmer is still needed to add new functions to the

library. Appendix B describes the procedure to follow.

As mentioned in the previous Chapters, at this stage
NOLMOD bond graphs are subject to the following
restrictions:

1. Storage and dissipation‘elements are one-ports.

2. Storage elements have integral causality.

3. Transformers and gyrators have constant moduli.
However, non-linear models often are more complicated. A
common situation is the need for a modulated transformer,
-MTF-, or modulated gyrator, -MGY-. Since the current
program can not handle such problems, further

implementations are necessary. The recommended sequence of

41

implementation is as follows :
1. Nonlinear transformers and gyrators.
2. Multi-port storages with integral causality.

-3. Multi-port storages with derivative causality.

LI ST OF REFERENCES

10.

11.

LIST OF REFERENCES

Frank H. Speckhart, Walter L. Green, A guide to using
CSMP - The Continuous System Modeling Program,
Prentice-Hall, 1976.

Continuous System Modeling Program III (CSMP III)
Program Reference Manual, SHl9-700l-3, IBM, 1975.

System/360 Continuous System Modeling Program User's
manual, GH20-0367-4, IBM, 1972.

DIFFEQ User's Manual, A. H. Case Center, College of
Engineering, Michigan State University, 1982.

James C. Bowers, John E. O'Reilly and Gary A. Show,
SUPER-SCEPTRE User's manual, University of South
Florida, Tampa, Florida, 1975.

Dean Karnopp and Ronald C. Rosenberg, System Dynamics:
A Unified Approach, John Wiley & Sons, 1975.

Ronald C. Rosenberg, ENPORT-5 USER'S MANUAL, A. H. Case
Center, College of Engineering, Michigan State
University, 1981.

H. M. Paynter, Analysis and Design of Engineering
Systems, M.I.T. Press, Cambridge, 1960.

S. H. Crandall, D. C. Karnoppr E. F. Kurtz, and
D. C. Pridmore-Brown, Dynamics of Mechanical and Electro
mechanical Systems, McGraw-Hill, 1968.

L. F. Godbout, Jr. and D. Jordan, On State Equation
Descriptions of Linear Differential Systems, ASME
Journal of Dynamic systems, Measurement, and Control,
Vol. 97, Dec. 1975.

R. C. Rosenberg, State-Space Formulation for Bond Graph

Models of Multiport Sysyems, ASME Journal of Dynamic
systems, Measurement, and Control, Vol. 93, March 1971.

42

43

12. Brice Carnahan, H. A. Luther, James O. Wilkes, Applied
Numerical Methods, John Wiley & Sons, 1969.

13. Herbert E. Merritt, Hydraulic Control Systems, John
Wiley & Sons, 1967.

14. Dean Karnopp, Using bond graphs in nonlinear simulation
problems, ASME Journal of Dynamic systems, Measurement,
and Control, Vol. 103, Dec. 1981.

APPENDI CES

APPENDIX A

HOW TO USE NOLMOD

APPENDIX A

HOW TO USE NOLMOD

A.l Introduction

NOLMOD is an interactive program for specifying nonlinear
constitutive equations that arise in the study of nonlinear dynamic

systems using bond graphs.

The program can be invoked from ENPORT by using the command
'FUNCTION'. Available commands include : LOAD, LIST, CHANGE, USER,
SYSTEM, HELP and RETURN. Each command can be abbreviated to the
smallest set of characters that maintains its uniqueness. After you
enter the command 'FUNCTION', NOLMOD will print out a function menu and
ask for a choice.

COMMANDS

FUNCTION

FUNCTION MENU :

ZBIB’T'EBXB'EQAPH mm W.

LIST : LIST OF FUNCTIONS TO BE SPECIFIED

CHANGE: CHANGE FUNCTION DEFINITION

USER : LOOK AT USER DEFINED FUNCTIONS

Al

A2

SYSTEM: LOOK AT SYSTEM FORMAT FUNCTIONS
HELP : HELP ABOUT FUNCTION
RETURN: RETURN TO MAIN MENU

ENTER OPTION : (L0, Ll, C, U, S, H, R ):

There are three additional commands that do not appear in the
function menu. These are INFO, DISP and EVAL. The set is used if you
wish to develop more understanding of the program operation. They are
used for helping the programmer doing further implementation. Usage of
the commands is described in next section in alphabetical order. The
normal processing sequence is : (DISP), LOAD, (INFO), LIST, (HELP),
CHANGE, USER or SYSTEM, (EVAL), RETURN. The commands shown in

parentheses are optional.

A.2 Commands and usage

I. CHANGE command

The CHANGE command lets you define or change the function type for
each equation by choosing functions from the function library. The
default function type is LINEAR. All the functions have been assigned
default parameters. The default value is shown in parentheses (). If
you accept the default value , just hit the return key. Otherwise,
enter the desired value. For a YES/NO question, if you type in the
wrong characters the system will treat it as NO. The name of the
function also can be abbreviated. You can terminate the operation of

CHANGE at the first three stages by entering 'X' in response to the

A3

question. Entering a zero '0' for function index and entering 'QUIT'
for a function type will also terminate the CHANGE process. Once you
have been asked for function parameters, you must continue the process.
No termination can be done. To leave quicker just accept the defaults.

The rules cited also hold for commands USER and SYSTEM. For example:

ENTER OPTION : (LO, LI, C, U, S, H, R ):
C

CHANGE ALL ? (NO):

ENTER FUNCTION INDEX :

I

FUNCTION I NOW DEFINED AS :
F I = LINEAR (P I)

REVERSE INPUT AND OUTPUT ? (NO):
Y

FUNCTION l NOW DEFINED AS :
P I - LINEAR (F I)

FUNCTION TYPE ?
SIN

P l = CI * SIN( C2 * F I + C3 ) + CA

VALUE OF Cl ? ( 1.0000 ):
SALUE OF C2 7 ( 1.0000 ):
3ALUE OF C3 ? ( 0.0000 ):
UALUE OF Ch ? ( 0.0000 ):
0.3I

VERIFY 7 (N0):
Y

P 1 = 2.00005+oo * SIN ( 3.00005+oo * F 1 + -h.OOOOE+OO )
+ 3.10005-01

CHANGE ANOTHER FUNCTION ? (YES):

Ah

YES

ENTER FUNCTION INDEX :
X

ENTER OPTION : (LO, LI, C, U, S, H, R ):

2. OISP command

This command turns on an extra print-switch that enables you to
view detailed information. It might help you to understand the internal

operation better.

3. EVAL command

This command is used to evaluate the functions that were defined
previously by the CHANGE operation. You can ask for a function index to
be evaluated and set the numerical value of the input variable of that
particular function. Note that all functions appear in system-formatted
form.

ENTER OPTION : (LO, LI, C, U, S, H, R ):
EVAL

GIVE FUNCTION NUMBER :
I

F l I 3.3333E-0l * ASIN( 5.0000E-OI * P I + -I.5500E-OI )
+ l.3333E+00

GIVE THE VALUE OF INPUT VARIABLE :
0.5

VALUE OF OUTPUT VARIABLE IS : I.3650E+00

EVALUATE ANOTHER FUNCTION ? (YES):
N

A5

A. FILE command

This is an ENPORT command used for saving data for later use. (See
ENPORT-5 User's Manual for a detailed description of ENPORT commands.)
Useful option for NOLMOD is:

FILE

GRAPH filename

The file name must not have more than four characters. The saved file
can be read back by NOLMOD's LOAD command or ENPORT's LOAD command,
which depends upon where you would like to restart it. However NOLMOD

only accepts a saved GRAPH file.
5. HELP command

The HELP command provides the list of different library functions
that you can use to define the functions. Detail descriptions will be
presented in the next section.

6. INFO command

The INFO command prints additional information concerning the

causality and power direction of the bond graph.

 

A6

7. LIST command

The LIST command prints out system-formatted equations at the
terminal. Input and output variables for C and | ports are expressed by
energy and co-energy variables: E, F, P and Q. Causal relationships
are incorporated. The input variable is on right of the equation, the
output variable is on the left side. Each function is labeled for
identification. The number after FCN is the function index. For

example:

ENTER OPTION : (LO, Li, C, U, S, H, R ):
LIST

INTEGRAL STORAGE ELEMENTS

E 3 = FCN 3 (Q 3)

Fl-FCNI(PI)

DISSIPATION ELEMENTS

E 5 = FCN 5 (F 5)

SOURCE ELEMENTS
E 2 8 FCN ( TIME )

F A = FCN ( TIME )

8. LOAD command

LOAD reads graph information from a previously saved file
containing a causal bond graph. Such a file can be created by the
ENPORT command FILE (see FILE command). If the user already has bond

graph data in ENPORT, the LOAD command should not be used. Every time

A7

the LOAD command is used common blocks will be initialized for accepting
new data, existing data will be destroyed. The User will be asked for a
file name after the LOAD command is given. If the file is loaded
successfully a message will shown on the screen. For example :

ENTER OPTION : (L0, Ll, C, U, S, H, R ):
LOAD

GIVE FILE NAME
SAMPLE

*** GRAPH FILE LOADING COMPLETED ***

9. RETURN command

The RETURN command put you back under ENPORT COMMANDS control.

l0. SYSTEM command

The SYSTEM command lists the system-formatted function for each

equation on the terminal. For example:

ENTER OPTION : (L0, Ll, C, U, S, H, R ):
S

LOOK AT ALL ? (NO):
ENTER FUNCTION INDEX :
l

FUNCTION l :

F l I 3.3333E-Ol * ASIN( 5.0000E-0l * P l + -l.5500E-Ol )
+ l.3333E+00

LOOK AT ANOTHER FUNCTION 2 (YES):
NO

AB

ll. USER command

The USER command lists the user-defined function for each equation

at the terminal. For example:

ENTER OPTION : (LO, LI, C, U, S, H, R ):
U

LOOK AT ALL ? (NO):
ENTER FUNCTION INDEX :
l

FUNCTION l

P I I 2.0000E+00 * SIN ( 3.0000E+00 * F I + -A.0000E+00 )
+ 3.lOOOE-0l

LOOK AT ANOTHER FUNCTION ? (YES):
N

A.3 Function library

The function library is the heart of the program. The user may
choose desired functions from the library. Diagrams of several
functions are shown in Figures A.l to Figure A.9. Meanings of the
parameter data are explained by the diagrams. Available functions are
described below, where Y is the output, X is the input, and Ci are
parameters of the function. A table arranged in alphabetical order is

shown in Table 2.

10.
II.

12.

13.
1A.

15.
I6.

17.

IS.

19.

20.

. CONST

. LINEAR

. POLY

. SIN

. COS

. ASIN

. ACOS

. DSIN

. DCOS

DASIN
DACOS
EXP
ALOG
A852
ABSQRT
RCPSQR

COULOM

BRKPT

SATUR

ADJSAT

Y I C0

Y I C0

Y I C0

Y I C1

Y I C1

Y I CI

Y I C1

Y I C1

Y I C1

Y I CI

Y I C1

Y I C1

Y I CI

Y I C1

Y I C1

Y I CO

Y I C1

Coulomb function.

Y I Y1 + C1*(X-X1),
Y I Y1 + C2*(X-XI),

Select one data point and define two slopes.

Y I Y1

3':

+

A

A9

CI*X

CIAX + c2*x**2 + C3*X**3 + Cthkkh

SIN(C2*X+C3) + CA
COS(C2*X+C3) + Ch
ASIN(CZIX+C3) + Ch
ACOS(C2*X+C3) + Ch
SIN(C2*X+C3) + CA
COS(C2*X+C3) + Ch
ASIN(C2*X+C3) + Ch
ACOS(C2*X+C3) + Ch
EXP(C2*X+C3) + Ch
ALOG(C2*X+C3) + Ch

X * ABS (X)

SIGN(X) A SQRT(ABS(X))

CI / ( 1+(X/C2)**2 )
SIGN(X)
(Fig. A.h)

X <IX1
X1< X

(Fig. A.l)

(Fig. A.2)

(Fig. A.3)

v - Yl + ((vI-Y2I/(xI-x2IIAIx-xi)I

Y I Y2

Saturation function.

..<

.<
IINI

Yl + Cl*(X-Xl)
Y] + ((Yl-Y2)/(Xl-X2))*(X-Xl),
Y2 + C2*(X-X2)

(Fig. A.6)

9

(Fig. A.S)

X <IXI
X1<X<X2
X2<IX

'X <IXI

x1<x<x2
X2<=X

21.

22.

PWLIN

INTPLO

A10

Select two data points and define two slopes. (Fig. A.7)
This function has great flexibility for approximating
various special functions. For example, it can be used
as an inverse of function 18 by defining very large
slopes at both ends. (Fig. A.8)

PIECEWISE LINEAR

Linear interpolation between data points. Maximum number
of data point is set to 20. User also defines the slopes
leading to the first point and from the last point.
INTERPOLATION POLYNOMIAL

Function 21 and 22 act as arbitrary function generators,
and both are based on Newton's interpolation polynomial
method. Function 21 is linear interpolation, and
function 22 is used as a high-degree interpolation
polynomial. The degree of interpolation can be 1, 2,
3,... up to 19. Usually a degree of two or three is
commonly used. End slopes are not supplied. The maximum
number of data points is 20. (Fig. A.9) The center part
of the interpolation results is more accurate than those
at the ends. Therefore it is helpful to extend the data

range to be wider than the range of interest.

A11

 

 

Figure A.l Function ABSZ

 

 

Figure A.2 Function ABSQRT

 

 

 

Figure A.3 Function RCPSQR

C1

 

—C1

 

Figure A.Q Function COULOM

(X2,Y2)

1 X

(X1,Y1)

 

Figure A.6 Function SATUR

c1 . C2
\\,/ \.
(x1.Y1) \,_,../
(X5aY5)X

 

Figure A.8 Function PWLIN

A12

Cl

Y C2

/ x

/(X1,y1)
/

 

Figure A.5 Function BRKPT

/
(W

C1

Y C2

M)X

 

Figure A.7 Function ADJSAT

(X1.

(X6,Y6)

./' '\. ./
X1)

 

 

Figure A.9 Function INTPLO

A13

Table 2. List of functions in alphabetic order

ABSQRT
ACOS

ADJSAT

ALOG
ASIN

BRKPT

CONST
COS
COULOM
DACOS
DASIN
DCOS

DSIN

DEFINITION
"""""" i'l'Ei'Ii'Iiééiii
Y = CI A SIGN(X) * SQRT(ABS(X))
Y I Cl * ACOS(C2*X+C3) + CA
Y . YI + C1*(X-Xl) , x <=x1
Y = YI + ((Yl-Y2)/(Xl-X2))*(X-Xl), X1<X<X2
Y . Y2 + C2*(X-X2) , x2<=x
Y . CI A ALOG(C2*X+C3) + CA
Y = CI * ASIN(C2*X+C3) + CA
Y = Y1 + Cl*(X-Xl), x <IXI
Y = Yl + C2*(X-X1), Xl< x
Y = C0
Y . CI A COS(C2*X+C3) + CA
Y - CI * SIGN(X)
Y = CI / ACOS(C2*X+C3) + CA
Y = CI / ASIN(C2*X+C3) + CA
Y . CI / COS(C2*X+C3) + CA
Y = CI / SIN(C2*X+C3) + CA
Y - CI * EXP(C2*X+C3) + Ch

EXP

INTPLO

LINEAR

POLY

PWLIN

INTERPOLATION POLYNOMIAL

Y I C0 + C1*X

Y - C0 + Cl*x + C2*X**2 + C3*X**3 + Ch*x**h

PIECEWISE LINEAR

RCPSQR

SATUR

SIN

Alh

Table 2. (cont'd)

DEFINITION

Y I C0 + C1 / (C2+X*A2)

Y - Yl , x <=x1
Y a YI + ((YI—Y2)/(XI-x2))*(x—XI). XI<x<x2
Y . Y2 , x2<=x

Y = CI A SIN(CzAx+C3) + Ch

APPENDIX B

IMPLEMENTATION OF A NEW FUNCTION

APPENDIX B

IMPLEMENTATION OF A NEW FUNTION

If a new function will be added to the library, changes must be
made to the following subroutines : FCNLIB, FCNLOOK, GETVAL and
HELPNEW. Source code for these subroutines can be found in Appendix C.
Also a new subroutine for the function is needed. Let's do an example
to show how to inplement a new function in the library. Suppose RCPSQR
is a new function to be added and there already are 21 functions in the
library. (Description of function RCPSQR can be found in Appendix A.)
Therefore the new one will be function number 22. This number will
serve as the identifier of the function type. The following

modifications should be made.

1. Subroutine FCNLIB

There are three important variables to be defined, namely, IFCNTP,
FCNSET and PARDF. Variable IFCNTP stores the function type number.
Therefore set IFCNTP(l)-22 for new function RCPSQR. FCNSET is an array
storing the name of the functions, therefore add the new name at the
22nd position of the list. (Data statement is in subroutine FCNLIB.)

PARDF is a 30-by-IO matrix which stores the default values of the

parameters of the functions. Hence, each function has ten storage
spaces available for its parameters. Unused spaces are filled with
zeros. The number of row stands for the function type number. RCPSQR
is the 22nd function in the library; therefore row 22 stores the
default parameters of the function RCPSQR. Set a new row of appropriate
default parameters for the new function. At last add a call statement
for the new subroutine:
ELSE IF(CMATCH(IFOUND).EQ.'RCPSQR')THEN
IFCNTP(I)=22

CALL FCNSU822(I)

2. Subroutine FCNLOOK

Just add two line of statement like these:
ELSE IF (IFCNTP(I).EQ.22) THEN

CALL FCNSU822(I)

3. Subroutine GETVAL

A few statements for calculating the output value (Y) are needed.
Add lines like these:
ELSE IF(NFCNTP.EQ.22)THEN
YOUTP-EQPAR(NOEQ,l)+EQPAR(NOEQ,2) / (1+(XINP/EQPAR(NOEQ,3))**2)
RETURN
where YOUTP is the output variable, XINP is the input variable and the
EQPAR's are the parameters of the function. The index NOEQ is the

equation number.

A. Subroutine HELPNEW

This subroutine prints library functions on the terminal when HELP
command is used. Add a line of print statement to describe the new
function.

PRINT *,'RCPSQR Y I C0 + C1 / (C2+X**2)'

5. New subroutine FCNSUBZZ

The subroutine performs the functions of reading parameters and
display of the equation. The very first thing that the subroutine
should do is to check whether the user asked for inversion or not. If
yes, and the function inversIon is not available, a warning message
should be displayed on the screen. The user should find another
alternative to fulfill his requirements. If the inversion is available,
some processes should be done which depend on the nature of the function
itself. The reader can refer to other function subroutines of the
program. Subroutine GENCON is used to accept input of constant
parameters C's. And subroutine GENPAR is used for reading a pair of
(X,Y) coordinates. Descriptions of the arguments for both subroutines

can be found in the source code in Appendix C.

APPENDIC C

SOURCE CODE FOR NOLMOD

APPENDIX C

SOURCE CODE FOR NOLMOD

Subroutines in this appendix appear in the following sequence :

I. MAINC 17. PARTAK 33. FCNSUBII
2. INIT 18. SUBWTR 3A. FCNSUBIZ
3. INFO I9. IOPVSR 35. FCNSUBI3
A. HELPNEW 20. GENPAR 36. FCNSUBIA
5. SWITZ 21. GENCON 37. FCNSUBIS
6. LOAD 22. COOREX 38. FCNSUBI6
7. CMREAD 23. FCNSUBI 39. FCNSUBI7
8. CMNCPR 2A. FCNSUBZ A0. FCNSU818
9. ERRMSG 25. FCNSUB3 AI. FCNSUB19
10. DEFINE 26. FCNSUBA A2. FCNSUBZO
11. FCNDEF 27. FCNSUBS A3. INTERP
12. EQRVS 28. FCNSUB6 AA. DTABLE
I3. SHOW 29. FCNSUB7 A5. FMIDV
1A. FCNWTR 30. FCNSU88 A6. FCNSUBZZ
15. FCNLIB 31. FCNSUB9 A7. NLEVAL
I6. FCNLOOK 32. FCNSUBIO A8. GETVAL

A9. COMFILE and COMAREA (common blocks)

C---MAINC ----------------------------------------------------

C
C
C

nonnnnnnnnnnnnnn

MAIN DRIVER
SUBROUTINE MAINC

WRITTEN BY: TSUN-YU CHOU, JUNE 1982, MICHIGAN STATE U.
LATEST UPDATE: FEB. 1983, TYC.

SOURCE PROGRAM FOR NONLINEAR SIMULATION-0F DYNAMIC SYSTEMS
BASED ON BOND GRAH THEORY. THIS PROGRAM CAN BE USED
INDEPENDENTLY OR FUNCTION AS A MODULE OF ENPORT PROGRAM.
(CALLED BY ENPORT BY THE COMAND ”FUNCTION")

OPTIONS CONTROL :
OPTl : (LOGIN, HELP )
(DEFINE, CHANGE, LOOK )
OPTZ : (L, U, s. N. x. Q)
OPT3 : (LINEAR, INTPLO )

COMMON BLOCKS -----

SINSERT COMFILE
SINSERT COMAREA

C ............................................................
DATA CMSETl/‘LOAD ','LIST ','USER ','SYSTEM','CHANGE',
+ 'HELP ','ZMNDO7','ZMNDOB','ZMNDOS',‘RETURN',
+ 'ZMNDll','ZMNDIZ','INFO ','DISP ','EVAL ',
+ 'ZMNDI6','ZMNDl7','ZMNDIB','ZMNDI9','ZMNDZO',
+ 'ZMNDZI','ZMNDZZ','ZMNDZ3',‘ZMNDZA','ZMNDZS',
+ 'ZMNDZ6','ZMNDZ7','ZMN028','ZMN029','ZMND30'/
c ............................................................
DATA ILUN,JLUN/l,l/
IPNT(I)=O
CALL INIT
CALL DEFINE
C ............................................................
Io OPTII'LOGIN '
CALL HELPNEw
lOO NRITE (JLUN,9000)
gooo FORMAT(/.'ENTER OPTION : (LD, Ll, C, U, s, H, R ): ')
CALL CMREAD
IF (NOCMN .EQ. 0) GO TO IO
CALL CMNCPR(CMSET1.30)
IF(IFOUND .EQ. 0) THEN
CALL ERRMSG(l)
GO TO IO
ENDIF
IF(|FOUND .GT. I) THEN
CALL ERRMSG(2)
GO TO 10
ENDIF
C
C ----- PROCESS -----

C

C

C---LOAD ----------
IF (CMATCH(IFOUND) .EQ. 'LOAD') THEN
CALL LOAD

c .

C---LIST ----------
ELSE IF (CMATCH(IFOUND) .EQ. ’LIST')THEN
OPTZI'L'
CALL DEFINE

C

C---USER ----------
ELSE IF (CMATCH(IFOUND) .EQ. 'USER')THEN
OPTZI'U'
CALL SHDN

C

C---SYSTEM ----------
ELSE IF (CMATCH(IFOUND) .EQ. 'SYSTEM')THEN
OPTZI'S'
CALL SHDN

C

C---CHANGE ----------
ELSE IF (CMATCH(IFOUND) .EQ. 'CHANGE')THEN
OPTZI'U'
CALL FCNDEF

C

C---HELP ----------
ELSE IF (CMATCH(IFOUND) .EQ. 'HELP') THEN
OPTI-‘HELP '
CALL HELPNEN

C

C——-RETURN ---------- -

ELSE IF (CMATCH(IFOUND .EQ. 'RETURN' ) THEN
C CALL EXIT

PRINT A,' '
RETURN

C

C---INFO ----------
ELSE IF(CMATCH(IFOUND) .EQ.'INFO')THEN
CALL INFO

C

C---DISP ----------
ELSE IF (CMATCH(IFOUND) .EQ. 'DISP') THEN
CALL szTz

C

C---EVAL ----------
ELSE IF (CMATCH(IFOUND) .EQ. 'EVAL')THEN
CALL NLEVAL

c

C---DTHERNISE ----------

C

ELSE
CALL ERRMSG(3)

C
ENDIF
C
c ............................................................
IF (NOCMN .LT. 100) GO TO 100
C
STOP
END
C
C---INIT -----------------------------------------------------
C
SUBROUTINE INIT
c
C INITIALIZATION
C
SINSERT COMAREA
C ............................................................
OPTZI'U'
C
DO 10 I=1,50
IFCNTP(I)=2
RVS(I)=0.

DO 10 JI1,10
EQPAR(I,J)IPARDF(2,J)
10 CONTINUE

RETURN
END
C
C---INFO -----------------------------------------------------
C
SUBROUTINE INFO
C
C PROVIDE ADDITIONAL INFORMATION ABOUT CAUSALITY AND POWER.

SINSERT COMFILE
SINSERT COMAREA

c ............................................................
DATA INFSET/‘CAUSAL','POWER ','HOME ','SWITOA','SWIT05',
+ 'SWIT06','SWIT07','SWIT08','SWIT09','SWITIO',
+ 'SWITII','SWIT12','SWITI3','SWITIA','SWIT15',
+ 'SWITI6','SWITI7','SWIT18','SWITI9','SWIT20'/
C ------------------------------------------------------------

lO WRITE(JLUN,20)

20 FDRMAT(/.'INFD : CAUSAL, POWER, HOME 2’/)
CALL CMREAD
IF(NOCMN.EQ.O)GOTO IO
CALL CMNCPR(INFSET,20)

IF(IFOUND.EQ.O)THEN
CALL ERRMSG(I)
GOTO 10

ELSE IF(IFOUND.GT.1)THEN
CALL ERRMSG(Z)

GOTO 10
ENDIF
C
C ----- CAUSAL ---------
C
IF(CMATCH(IFOUND) .EQ.'CAUSAL') THEN
C
NBD2=NDDA2
DO 50 K=l,NBDZ
IF(ICMx(K).CE.5) NBIMX(K)=-NBIMX(K)
50 CONTINUE
WRITE (JLUN,SOlO)
NPI=I
DO 52 I: I.NEL
NP2=NPTR(I+l)-l
WRITE (JLUN,5020) |,IELNAM(I),(NBIMX(J),J=NP1,NP2)
52 NPl=NP2+1
C
00 5A K=1.NBDZ
5A NBIMX(K)=IABS(NBIMX(K))
GOTO 10
C
50I0 FORMAT(lX/' THE BOND GRAPH STRUCTURE WITH CAUSALITY IS'
+ ,/,3x,'(+ = STROKE END, - = NON STROKE END)'/)
5020 FORMAT(2X,13,2x,A2,2X,A(15IA/9X))
C
C ----- POWER ------------
C
ELSE IF(CMATCH(IFOUND).EQ.'POWER')THEN
C
WRITE (JLUN,lOOO)
N2=0
101 NI=N2+I
N2=NBD
IF ((N2-Nl).GT.9) N2INI+9
wRITE(JLUN,I0I0) (N,N=N1,N2)
WRITE(JLUN,1020) (IELNAM(IBMX(N,1)),NIN1,N2)
WRITE(JLUN.1030) (IELNAM(IDMX(N,2)).N=NI,N2)
IF (N2.LT.NBD) GOTO 101
GOTO I0
C
I000 FORMAT(1X/' THE POWER DIRECTIONS‘)
lOlO FORMAT(lX/' BOND'.2x.I2.9(3x,I2))
I020 FORMAT(' FROM',2X,AA,9(1X,AA))
1030 FORMAT(' TO',2X,AA,9(IX,AA))
C
C ----- HOME ----------
C
ELSE IF(CMATCH(IFOUND).EQ.'HOME')THEN
RETURN
C
C

C ----- OTHERHISE ---------

C

ELSE

CALL ERRMSG(3)
GOTO 10

C

ENDIF
C

RETURN

END
C
C---HELPNEW ---------------------------------------------------
C

SUBROUTINE HELPNEw
C
C OPTl ( HELP, LOGIN )

SINSERT COMFILE
SINSERT COMAREA
C ............................................................

C

C

PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT

IF(OPT1.EQ.'HELP ')GOTO 20
:‘<,' '
A,'FUNCTION MENU :'
S':,' --------------- '
*,I I
*,'LOAD : LOAD GRAPH FROM FILE'
*,'LIST : LIST OF FUNCTIONS TO BE SPECIFIED'
*,'CHANGE: CHANGE FUNCTION DEFINITION'
*,.USER : LOOK AT USER DEFINED FUNCTIONS'
*,'SYSTEM: LOOK AT SYSTEM FORMAT FUNCTIONS'
*,'HELP : HELP ABOUT FUNCTION'
*,'RETURN: RETURN TO MAIN MENU'
:‘t,' '

PRINT

IF(OPT1.EQ.'LOGIN ‘)RETURN

PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT

{"1 I

A,'NHEN USING COMMANDS “CHANGE”, ”USER”, ”SYSTEM“,'
A,'Y0U ALWAYS CAN TERMINATE THE PROCESS BY ENTER "x"'
A,'F0R THE FOLLOWING QUESTIONS :'

A,' (I) CHANCE ALL? DR LOOK AT ALL? '

A,' (2) ENTER FUNCTION INDEX? '

*.' (3) REVERSE INPUT AND OUTPUT?‘

A,'ENTER "0" FOR FUNCTION INDEx AND ENTER'

A,'"QUIT" FOR FUNCTION TYPE ALSO TERMINATE'

A,'THE PROCESS AND PUT YOU BACK TO MAIN MENU.‘

*’1 I

A,'HIT RETURN KEY TO CONTINUE.‘

A,'(IF YOU ARE USING TEKTRDNIx TERMINAL.‘

A,' PLEASE PACE THE SCREEN FIRST.)'

CALL CMREAD
PRINT *,'FUNCTION LIBRARY :'
PRINT *,' '

PRINT

' PRINT

PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT

PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT

PRINT *

PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT
PRINT

*,'TYPE

*,' ----------

*,'CONST
,0"! I
*,'LINEAR
*’I I
*,'POLY
*,l I
*,.SIN
,fi’l I
*,'ASIN
*,l I
*,'COS
,o": I
*,'ACOS
*,I I
*,'DSIN
*,I I
*,'DASIN
,9"I I
*,'DCOS
:‘F,’ '
*,'DACOS
*,I I
A,'EXP
*’I I
*,'ALOG
*’I I

DEFINITION'

Y = C0 + CIAx'

Y = C0 + CIAx + C2AXAA2 + C3AXAA3 + CAAXAAA'
Y = CI A SIN(C2Ax+C3) + CA'

Y = CI A ASIN(C2Ax+C3) + CA'

Y = CI A COS(C2Ax+C3) + CA'

Y a CI A ACOS(C2Ax+C3) + CA'

Y = CI SIN(C2Ax+C3) + CA'
Y = CI ASIN(C2Ax+C3) + CA'

Y I C1 COS(C2*X+C3) + CA'

\\\\

Y I C1 ACOS(C2*X+C3) + CA'
Y I C1 * EXP(C2*X+C3) + CA'

Y I C1 * ALOG(C2*X+C3) + CA'

*,'*** MORE, PAGE AND HIT (RETURN) KEY.‘
CALL CMREAD

*,I I
*,'TYPE

A,' ----------

*,'ABSZ
*’I I
*,'ABSQRT
3","
A,'RCPSQR
*’I I
*.'COULOM
*’I I
*,'BRKPT
*,'

*,'SATUR

*I
*I
*II
A 'ADJSAT
I
I
I
I

O
9
9
9

:5-

2':

at

*, PWLIN

DEFINITION'

Y I CI * X * ABS(X)'
Y I CI * SIGN(X) * SQRT(ABS(X))'
Y I C0 + C1 / (C2+X**2)'

Y = CI A SIGN(X)'

Y . Y1 + CIA(x-x1), x $8X1'

Y = Y1 + C2A(x-XI), Xl< x'

Y = YI , x <-Xl'
Y = Yl + ((Yl-YZI/(Xl-X2))*(X-Xl), Xl<x<X2'
Y . Y2 , x2<-x'

Y - Y1 + CIA(x-XI) , x <=X1'
Y . YI + ((YI-Y2)/(XI-x2))A(x-XI), Xl<X<X2'
Y = Y2'+ C2A(x-x2) , x2<=x'

PIECEWISE LINEAR'

PRINT A,' '

PRINT A,'INTPLO INTERPOLATION POLYNOMIAL'
C

RETURN

END
C
C---SWIT2 ----------------------------------------------------
C

SUBROUTINE SNIT2
C
C SWITCH ROUTINE-- FOR ADDITION INFORMATION DISPLAY.
C ONLY I SWITCH IS USED (IPNT(I)) RIGHT NOW. UP TO 10
C SWITCHS CAN BE USED FOR FUTURE EXPANSION.
C lION, OIOFF.

$INSERT COMFILE
$INSERT COMAREA

c ............................................................
ITRANSIIPNT(I)
IF(ITRANS.EQ.O) THEN
IPNT(1)I1
ELSE IF(ITRANS.EQ.1) THEN
IPNT(1)=0
ENDIF
C
RETURN
END
C
C---LOAD -----------------------------------------------------
C
SUBROUTINE LOAD
C
C ----- WRITTEN BY : TSUN-YU CHOU, JULY 1982
C
C ----- LOAD DATA FROM SAVED FILE FOR FURTHER PROCESS.(BY SUB DEFINE)
C DATA FILE CONTAIN NECESSARY GRAPH DATA SAVED BY
C ENPORT COMMAND :
C FILE
C GRAPH filaname
C
LOGICAL EXIST, OPENED
CHARACTER*BO FILNAM
C

$INSERT COMFILE
$INSERT COMAREA
C ............................................................

C
c .....
C

c .....

INITIALIZE THE PARAMETER LIST.

CALL INIT
CALL CMNZER

READ FILE NAME.

10

WRITE(JLUN,lOOO)

CALL CMREAD
IF(NOCMN.EQ.O)GOTOIO
IF(CMN.EQ.':') RETURN
FILNAMICMN

CHECK IF FILE EXIST.A

INQUIRE(FILE-FlLNAM,ExIST-EXIST)
|F(EXIST .EQV..FALSE.)THEN
WRITE(JLUN,IOOI)

RETURN

FILE EXIST, SEARCH FOR UNOPENED UNIT.

ELSE
NUNITIS

INQUIRE(UNITINUNIT,OPENEDIOPENED)

IF(OPENED.EQV..TRUE.)THEN
NUNITINUNIT+I

GOTO 30

ELSE
OPEN(NUNIT,FILEIFILNAM)
ENDIF

ENDIF
HEADING INFO.

READ(NUNIT,1010,ERRI100,ENDI100)
READ(NUN|T,1010,ERRI100,ENDI100)
READ(NUNIT,IOAO,ERRI100,ENDI100)

GRAPH DATA BLOCK.

READ(NUNIT,I0I0.ERR-I00,END-IOO)
READ(NUNIT,10AO,ERRI100.ENDIIOO)
READ(NUNIT,I0I0.ERR-I00,END=IOO)
READ(NUNIT,1020,ERRIIOO,ENDIIOO)
I2-2ANDD

READ(NUNIT,I020,ERR=I00,END-I00)
READ(NUNIT.I020.ERR=I00.END=I00)
IZINEL+I

READ(NUNIT,I020.ERR=IOO,END-I00)

PARAMETERS

I2-NEL+I

READ(NUNIT,IOIo,ERR=I00,END=IOO)
READ(NUNIT,1020,ERR=lOO,END=lOO)
READ(NUNIT,I050,ERR=IOO.END=IOO)

NAMEF
(HEAD(I),III,20)
ILUN.JLUN

KNORD

NEL,NBD
(IELNAM(I),I=1,NEL)
(IELLST(I),III,NEL)

(NBIMX(I),III,IZ)
((IBMX(I,J),JI1,2),III,NBD)

(NPTR(I),I-I,I2)

KNDRD
(IPTR(I),II1,I2),NPAR.IPFLG
(PAR(I),II1,NPAR)

c .....
C

100

CAUSALITY

READ(NUNIT,1010,ERRI100,ENDI100) KWORD

12I2*NBD

READ(NUNIT,I020,ERR=IOO,END-I00) (ICMX(|),III,12)
READ(NUNIT,I020.ERR=I00,END-IOO) (IBEQ(I),I=1,NBD)
READ(NUNIT,I020,ERR-I00,END=IOO) (IDT(I).I=I,NBD)
READ(NUNIT,IOAO,ERR=100,ENDIIOO) NBEX,NFI,NFD,NFL,NFS,NBIN

CLOSE(NUNIT,STATUSI’KEEP')

ECHO OF READING

IF (IPNT(I).EQ.O) GOTO 50

WRITE(JLUN,1010)
WRITE(JLUN,1010)
WRITE(JLUN,10AO)

WRITE(JLUN,2010)
WRITE(JLUN,10AO)
WRITE(JLUN,1010)
WRITE(JLUN,1020)
I2=2ANDD

WRITE(JLUN,1020)
WRITE(JLUN,1020)
I2=NEL+I

WRITE(JLUN,1020)

PARAMETERS

IZINEL+I

WRITE(JLUN,2020)
WRITE(JLUN,1020)
WRITE(JLUN,1050)

CAUSALITY
wRITE(JLUN.2030)
I2=2ANDD
WRITE(JLUN,1020)
WRITE(JLUN,1020)
WRITE(JLUN,1020)
WRITE(JLUN,10AO)
wRITE(JLUN.900)
CALL DEFINE

RETURN

NAMEF
(HEAD(I),III,20)
ILUN.JLUN

NEL.NBD
(IELNAM(I),|II,NEL)
(IELLST(I),I=I,NEL)

(NBIMX(|),II1,12)
((IBMX(I,J),JI1,2),III,NBD)

(NPTR(I),III,|2)

(IPTR(I),III,12),NPAR,IPFLG
(PAR(I),III,NPAR)

(ICMX(|),I=1,I2)
(IBEQ(I),I=I,NBD)
(IBT(I).I=1.NBD)
NBEX,NFI,NFD,NFL,NFS,NBIN

WRITE(JLUN.990) FILNAM
CLOSE(NUNIT,STATUSI'KEEP')

RETURN
C
C ----- FORMAT ----------
C

900 FORMAT(/.'*** GRAPH FILE LOADING COMPLETED ***')

990 FORMAT(/,'OOPS ] YOU ARE LOADING THE WRONG FILE : ',A80)
1000 FORMAT(/,'GIVE FILE NAME :')

1001 FORMAT(/,'FILE NOT EXIST')

lOlO FORMAT(ZOAA)

I020 FORMAT(20I3)

lOAO FORMAT(6IA)

I050 FORMAT(5E12.5)

20I0 FORMAT(‘GRAPH')

2020 FORMAT(‘PARAMETERS')
2030 FORMAT('CAUSAL')

C
END
C
C---CMREAD ---------------------------------------------------
C
SUBROUTINE CMREAD
C
C ----- WRITTEN BY : TSUN-YU CHOU, JUNE 1982
C
C ----- FREE FIELD COMMAND READER
C

$INSERT COMFILE
$INSERT COMAREA

C
INTEGER BLANK
LOGICAL FOUND
DATA BLANK/160/
c ............................................................

100 CONTINUE
READ(JLUN,9000) CINPUT
FPOCMNIO
LPOCMNIO
FOUNDI .FALSE.

DO 500 I-I,80
IF( ICHAR( CINPUT(I:I) ) .EQ. BLANK) THEN
IF(FPOCMN .EQ. 0) THEN
GO TO 500
ELSE IF(FPOCMN .NE. 0) THEN
LPOCMNIl-l
NOCMNII
GO TO 800
ENDIF
ELSE
IF(FOUND .EQV. .TRUE.) GO TO 500
FPOCMN=I
FOUND= .TRUE.

ENDIF
500 CONTINUE

IF((FPOCMN .EQ. O ) .AND. (LPOCMN .EQ. O) ) THEN
NOCMNIO
RETURN

ELSE IF((FPOCMN .NE. 0) .AND. (LPOCMN .EQ. 0) ) THEN
WRITE(JLUN,9OOI)
NOCMNIO
GO TO 900

ENDIF

800 CONTINUE
CMNICINPUT(FPOCMN:LPOCMN)

c ----- FORMAT
9000 FORMAT(ABO)
9001 FORMAT(//.'YOUR COMNAND EXCEED 80 COLUMN ')

900 RETURN
END
C
C---CMNCPR --------------------------------------------------
C
SUBROUTINE CMNCPR(CMNSET,NOCMAP)
C
C---WRITTEN BY : TSUN-YU CHOU. JUNE 1982
C
C THIS SUBROUTINE MAPPING THE INPUT COMMAND TO THE
C COMMAND LIST, COMPARING AND ACCEPTING ABBRIVATED
C COMMAND.
C
CHARACTER*6 CMNSET(100)
C

$INSERT COMFILE
$INSERT COMAREA

c ............................................................
NOWORDILPOCMN-FPOCMN+I
IFOUNDIO

C
00 200 II1,NOCMAP

C

D0 100 JI1,NOWORD
IF ( CMNSET(I)(J:J) .NE. CMN(J:J) ) GO TO 200
100 CONTINUE

C
IFOUNDIIFOUND+1
CMATCH(IFOUND)ICMNSET(I)
C
200 CONTINUE
C

RETURN
END

C
C ----- ERRMSG -------------------------------------------------
C
SUBROUTINE ERRMSG(IE)
C
C THIS SUBROUTINE PRINT OUT (ERROR) MESSAGE IN THE TERMINAL.
C

$INSERT COMFILE
$INSERT COMAREA
c ............................................................

1001
1002
1003
100A
1005
1006

1007
1008
1009

+

GOTO(10,20,30,A0,50,60,70,80,90,100),IE

WRITE(JLUN,1001)CMN
RETURN

WRITE(JLUN,1002)(CMATCH(K),KI1,IFOUND)
RETURN

WRITE(JLUN,1003)CMN
RETURN

WRITE(JLUN.100A)CMN
RETURN

WRITE(JLUN,1005)
RETURN

wRITE(JLUN.1006)
RETURN.

WRITE(JLUN,1007)
RETURN

WRITE(JLUN,1008)
RETURN

WRITE(JLUN,1009)
RETURN

RETURN
FORMAT ------------

FORMAT(/,' ERROR] NO SUCH COMMAND : ',A80)
FORMAT(/,' ERROR] COMMANDS MIXED UP',/,(A6,'?'))
FORMAT(/.' SORRY] THIS COMMAND CANNOT USE RIGHT NOW : ',A80)
FORMAT(/,' ERROR] EXPECTING NUMBER, YOUR INPUT IS : ',A80/)
FORMAT(/.' ERROR] BAD PARAMETER. ',/)
FORMAT(/.' YOU MAY NOT REVERSE CAUSALITY FOR THIS FUNCTION',

/,' EITHER CHANGE THE FUNCTION OR ACCEPT THE CAUSALITY.')
FORMAT(/,'VERIFY ? (NO):')
FORMAT(/,' ERROR] INVERSE OF ZERO SLOPE WILL BE INFINITY.',/)
FORMAT(/,' ERROR] SLOPE CAN NOT BE INFINITY.'/)

C
END
C
C---DEF|NE ---------------------------------------------------
C
SUBROUTINE DEFINE
C
C ----- WRITTEN BY : TSUN-YU CHOU. JULY I982
C

C---PURPOSE: DEFINE THE CONSTITUTIVE EQUATIONS OF BG.
C

C ----- VARIABLES:

C IBT(I) --- CLASSIFY TYPES OF BOND
C POSITIVE-- E INTO FIELD ELEMENT
C NEGATIVE-- F INTO FIELD ELEMENT
C

C TYPES OF BONDS IBT

C EXTERNAL

C INDEPENDENT I

C DEPENDENT 2

c R 3

C SOURCE A

C INTERNAL 5

C

C IEELST --— TYPES OF ELEMENTS.

C TYPES OF ELEMENTS

C C I

C I 2

c R 3

C SE A

C SF 5

C 0 6

C I 7

C TF 8

C GY 9

C MACRO 10

C

C

C NEQTP --- TYPES OF EQUATIONS.

C TYPE EQUATION

c ............................
c I EI - FCN ( Q1 )
C

C 2 Fl = FCN ( PI )
C

C 3 QI - FCN ( EI )
C

c A PI - FCN ( FI )
C

C 5 EI = FCN ( Fl )
C

C 6 Fl = FCN ( EI )
C

C 15

----- COMMON BLOCK

C 7 E1 = FCN (TIME)

C

C 8 Fl - FCN (TIME)

C

c 12 EI = FCN ( E2 )

C F2 = FCN ( Fl )

C

C I3 FI . FCN ( F2 )

C E2 = FCN ( EI )

C

C IA EI . FCN ( F2 )

C E2 = FCN ( FI )

C

C 15 F1 = FCN ( E2 )

C F2 = FCN ( El )

C

C NFCN --- STORE THE FUNCTION IDENTIFIER.
C NOIO --- TOTAL NUMBER OF I/D VARIABLES.
C IEQZS(NEQS.I) --- OUTPUT PORT IDENTIFIER.
C IEQ2$(NEQS,2) --- INPUT PORT IDENTIFIER.
C CHIO(NEQS,l/3) --- OUTPUT VARIABLE.

C CHIO(NEQS,2/A) --- INPUT VARIABLE.

C

C

C

DIMENSION IGNORE(50)
C
$INSERT COMFILE
$INSERT COMAREA
c ............................................................
C ----- CHECKING FOR HEADING.
C
OPTII'DEFINE'
IF(IPNT(I).EQ.I) OPTZI'L'
IBTc-O
IBTIIO
IBTRIO
IBTSIO
DO 5 I-I,NBD
IBTABS-IABS(IBT(I))
IF(IBTABS.EQ.I) IBTC-IBTC+I
IF(IBTABS.EQ.2) IBTI-IBTI+I
IF(IBTABS.EQ.3) IBTR-IBTR+I
IF(IBTABS.EQ.A) IBTS=IBTS+I
CONTINUE

----- CHECKING FOR ONE-PORT C. I. R ELEMENTS.
DO 20 II1,NBD
IGNORE(I)I1

20 CONTINUE

NPIII

DO 30 II1,NEL
NP2INPTR(I+1)-I
IF((IELLST(I).LE.3).AND.(NP1.NE.NP2))GOTO 25
GOTO 30

25 DO 27 KINPI,NP2
IGIIABS(NBIMX(K))

27 IGNORE(IG)I0
30 NPIINP2+I
C
C ----- PRINT OUT HEADING AND SYMBOLIC EQUATIONS FOR
C ----- ONE-PORT C, I. R ELEMENTS AND SOURCE ELEMENTS SE AND SF.
C
NEQSIO
C
00 200 II1,A

IF(OPT2.EQ.'L')THEN
IF((I .EQ. l).AND.(IBTC.GT.O)
IF((I .EQ. 2).AND.(IBTI.GT.0)
IF((I .EQ. 3).AND.(IBTR.GT.0)
IF((I .EQ. A).AND.(IBTS.GT.O)

ENDIF

ICAUS=-IAI

) WRITE(JLUN,9001)
) WRITE(JLUN,9002)
) WRITE(JLUN,9003)
) WRITE(JLUN,900A)

50 DO I00 J=1.NBD
IF(IBT(J) .NE. ICAUS) GOTOIOO
IF((IBT(J).EQ. l).AND.(|GNORE(J).EQ.l)) GOTO 65
IF((IBT(J).EQ.-l).AND.(IGNORE(J).EQ.1)) GOTO 60
IF((IBT(J).EQ. 2).AND.(IGNORE(J).EQ.1)) GOTO 70
IF((IBT(J).EQ.-2).AND.(IGNORE(J).EQ.l)) GOTO 75
IF((IBT(J).EQ. 3).AND.(IGNORE(J).EQ.1)) GOTO 85
IF((IBT(J).EQ.-3).AND.(IGNORE(J).EQ.l)) GOTO 80
IF((IBT(J).EQ. A).AND.(IGNORE(J).EQ.I)) GOTO 95
IF((IBT(J).EQ.-A).AND.(IGNORE(J).EQ.1)) GOTO 90
C
GOTO I00
C
C ----- INTEGRATE C-ELEMENT
60 NEQS=NEQS+I
NFCN(NEQS)-J
NEQTP(NEQS)-I
IEQZS(NEQS,1)=J
IEQ2$(NEQS,2)=J
N0I0(NEQS)=2
CHIO(NEQS,1)='E'
CHIO(NEQS,2)-'Q'
IF(OPT2.EQ.'L')THEN
CALL FCNWTR(NEQS)
ENDIF
GOTO I00
C
C ----- INTEGRATE I-ELEMENT
65 NEQS=NEQS+I
NFCN(NEQS)=J

NEQTP(NEQS)I2
IEQ2$(NEQS,I)=J
IEQZS(NEQS,2)IJ
NOIO(NEQS)=2
CHIO(NEQS,l)-'F'
CHIO(NEQS,2)='P'
IF(OPT2.EQ.’L')THEN
CALL FCNWTR(NEQS)
ENDIF
GOTO I00
C
C ----- DERIVATIVE C-ELEMENT
70 NEQS=NEQS+I
NFCN(NEQS)=J
NEQTP(NEQS)=3
IEQZS(NEQS.1)IJ
IEQZS(NEQS,2)IJ
NOIO(NEQS)=2
CHIO(NEQS,l)-'Q'
CHIO(NEQS,2)-'E'
IF(OPT2.EQ.'L')THEN
CALL FCNWTR(NEQS)
ENDIF
GOTO 100
C
C ----- DERIVATIVE I-ELEMENT
75 NEQS=NEQS+I
NFCN(NEQS)IJ
NEQTP(NEQS)-A
IEQ2$(NEQS,I)=J
IEQZS(NEQS,2)=J
NOIO(NEQS)=2
CHIO(NEQS,I)='P'
CHIO(NEQS,2)='F'
IF(OPT2.EQ.'L')THEN
CALL FCNWTR(NEQS)
ENDIF
GOTO 100
C
c ----- FLOW INPUT To R ELEMENT
80 NEQS=NEQS+I
NFCN(NEQS)IJ
NEQTP(NEQS)-5
IEQZS(NEQS,I)=J
IEQZS(NEQS.2)IJ
NOIO(NEQS)-2
CHIO(NEQS,l)-'E'
CHIO(NEQS,2)-'F'
IF(OPT2.EQ.'L')THEN
CALL FCNWTR(NEQS)
ENDIF
GOTO lOO

C ----- EFFORT INPUT TO R ELEMENT

85 NEQS=NEQS+I
NFCN(NEQS)=J
NEQTP(NEQS)-6

IEQ2$(NEQS,I)-J
IEQZS(NEQS,2)-J
NOIO(NEQS)IZ
CHIO(NEQS.I)='F'
CHIO(NEQS,2)I'E'
IF(OPT2.EQ.'L')THEN
CALL FCNWTR(NEQS)

ENDIF
GOTO I00
C
C ----- SOURCE ELEMENT SE
90 IF(OPT2.EQ.'L')THEN
WRITE(JLUN,9OIO) J
ENDIF
GOTO I00
C
C ----- SOURCE ELEMENT SF
95 IF(OPT2.EQ.'L')THEN
WRITE(JLUN.9020) J
ENDIF

C
100 CONTINUE
C

IF(ICAUS .GT. 0)GOTO 200

ICAUSII
GOTO 50
200 CONTINUE
C
C ----- SYMBOLIC EQUATIONS FOR MULTl-PORT ELEMENTS
C
NPlII
DO 700 I-I,NEL
NP2-NPTR(I+I)-I
IF(NPI .EQ. NP2)GOT0 700
NP-NP2—NPI+I
GOTO(610,620,630,700,700,700,700,700,700,700),IELLST(I)
C
C ----- MULTI-PORT C-FIELD
C

610 IF(OPT2.EQ.'L') WRITE(JLUN,8100)NP
DO 611 J-1,NP
IF(OPT2.EQ.'L')THEN
WRITE(JLUN,8110) NBIMX(J),NBIMX(J),(NBIMX(K),KINP1,NP2)
ENDIF
NEQS-NEQS+1
NFCN(NEQS)INBIMX(J)
NEQTP(NEQS)-9
611 CONTINUE
GOTO 700

C
C ----- MULTI-PORT I-FIELD
C
620 IF(OPT2.EQ.'L') WRITE(JLUN,8200)NP
DO 621 J-1,NP
IF(OPT2.EQ.'L')THEN
WRITE(JLUN,82IO) NBIMX(J),NBIMX(J),(NBIMX(K),KINP1,NP2)
ENDIF
NEQSINEQ+1
NFCN(NEQS)-NBIMX(J)
NEQTP(NEQS)IIO
621 CONTINUE
GOTO 700
C
C ----- MULTI-PORT R-FIELD
C
630 IF(OPT2.EQ.'L') WRITE(JLUN,8300)NP
DO 631 JIl,NP
IF(OPT2.EQ.'L')THEN
WRITE(JLUN,83IO) NBIMX(J),NBIMX(J),(NBIMX(K),K=NP1,NP2)
ENDIF
NEQS=NEQS+I
NFCN(NEQS)-NBIMX(J)
NEQTP(NEQS)=II
631 CONTINUE
C
700 NPI-NP2+I
C
C ----- SYMBOLIC EQUATIONS FOR TF AND GY.
C
ISLT-B
690 NPI=I
DO 800 I-I,NEL
NP2=NPTR(I+I)-I
IF(NPI .EQ. NP2)GOT0 800
NP-NP2-NPI+I
IF((IELLST(I) .EQ. 8) .AND. (ISLT .EQ. 8))GOTO710
IF((IELLST(I) .EQ. 9) .AND. (ISLT .EQ. 9))GOTO720
GOTO 800
C
C ----- TRANSFORMER
C
710 CONTINUE
IF(ICMX(NPI) .LE. 3)GOT0 7II
C

IF(OPT2.EQ.'L') WRITE(JLUN,8800)
NEQS=NEQS+I
NFCN(NEQS)-NBIMX(NPI)
NEQTP(NEQS)-I2
IEQZS(NEQS,I)INBIMX(NPI)
IEQZS(NEQS,2)INBIMX(NP2)

NOIO(NEQS)IA
CHIO(NEQS,1)='E'

CHIO(NEQS,2)I'E'
CHIO(NEQS,3)I'F'
CHIO(NEQS,A)I'F'
IF(OPT2.EQ.'L')THEN
CALL FCNWTR(NEQS)
ENDIF

GOTO 800
C

7II IF(OPT2.EQ.'L') WRITE(JLUN,8800)
NEQS=NEQS+I
NFCN(NEQS)INBIMX(NPl)
NEQTP(NEQS)=I3
IEQZS(NEQS,1)INBIMX(NP1)
IEQZS(NEQS,2)INBIMX(NP2)
NOIO(NEQS)IA
CHIO(NEQS,1)I'F'
CHIO(NEQS,2)-'F'
CHIO(NEQS,3)I'E'
CHIO(NEQS,A)I'E'

IF(OPT2.EQ.'L')THEN
CALL FCNWTR(NEQS)
ENDIF
GOTO 800

C ----- GYRATOR

720 CONTINUE
IF(ICMX(NPI) .LE. 3)GOTO 721

C
IF(OPT2.EQ.'L') WRITE(JLUN.8900)
NEQS=NEQS+I
NFCN(NEQS)-NBIMX(NPI)
NEQTP(NEQS)IIA
IEQ2$(NEQS,I)-NBIMX(NPI)
IEQ2$(NEQS,2)=NBIHX(NP2)
NOIO(NEQS)=A
CHIO(NEQS,1)I'E'
CHIO(NEQS,2)=‘F'
CHIO(NEQS,3)I'E'
CHIO(NEQS,A)-'F'
IF(OPT2.EQ.'L')THEN
CALL FCNWTR(NEQS)
ENDIF
GOTO 800
C
72I IF(OPT2.EQ.'L') WRITE(JLUN.8900)
NEQS=NEQS+I
NFCN(NEQS)-NBIMX(NPI)
NEQTP(NEQS)-I5

IEQZS(NEQS,1)=NBIMX(NPI)

lEQ2$(NEQS,2)INBIMX(NP2)
NOIO(NEQS)IA
CHIO(NEQS,1)='F'

CHIO(NEQS,2)-'E'
CHIO(NEQS,3)I'F'
CHIO(NEQS,A)-'E'
IF(OPT2.EQ.'L')THEN
CALL FCNWTR(NEQS)
ENDIF
GOTO 800

800 NPI-NP2+I
IF(ISLT .EQ. 9)GOTO 810

ISLTI9
GOTO 690
C
810 CONTINUE
C
IF (IPNT(I).EQ.0) GOTO 1001
C

WRITE(JLUN,9IOO)

DD I000 I=I,NEQS

WRITE(JLUN,9200)I,NFCN(I),NEQTP(I)
lOOO CONTINUE

C
1001 OPT2I'U'
RETURN
C
C--- FORMATS
C
C ----- CLASSIFIED HEADING
C

9OOI F0RMAT(//.25HINTEGRAL STORAGE ELEMENTS)
9002 FORMAT(//.27HDERIVATIVE STORAGE ELEMENTS)
9003 FORMAT(//.20HDISSIPATION ELEMENTS)

900A FORMAT(//.I5HSOURCE ELEMENTS)

C ----- SE AND SF

C

90I0 FORMAT(/,5X,'E',I2,' a FCN ( TIME )')

9020 FORMAT(/,5X,'F',12,' = FCN ( TIME )')

C

C ----- MULTI-PORT ELEMENTS

C

8IOO FORMAT(//.12,1AH-PORT C—FIELD )

8110 FORMAT(/,5X,'E',12,' = FCN',12,' (Q',12,(', Q',12),')' )
C

8200 FORMAT(//.12.1AH-PORT I-FIELD )

8210 FORMAT(/,5X,'F',|2,' - FCN',12,' (P',I2,(', P',12),‘)' )
C

8300 FORMAT(//.I2.IAH-P0RT R-FIELD )

83I0 FORMAT(/,5X,'E',|2,' = FCN',I2,' (F',|2,(', F',|2),')' )
C

C ----- TF AND GY
C
8800 FORMAT(//.IIHTRANSFDRMER)

8900 FORMAT(//.7HGYRATOR)

9I00 FORMAT(/5x,3HNO.,5X,3HFCN,5X,AHTYPE)
9200 FORMAT(5X,I2,6X,I2,7X,I2)

C
END
C
C
C-'-FCNDEF ---------------------------------------------------
c.
SUBROUTINE FCNDEF
C
C ----- WRITTEN BY : TSUN-YU CHOU. JULY 1982
C
C OPT2 I N ---> WHEN FUNCTION INVERSE NOT ABAILABLE.
C OPT2 I Q ---> USER ENTER QUIT FOR FUNCTION TYPE.
C 0PT2 I X ---> USER ENTER X TO TERMINATE THE OPERATION.

$INSERT COMFILE

$INSERT COMAREA

c ............................................................
OPTII'CHANGE'

I25 WRITE(JLUN.I26)
I26 FORMAT(/,'CHANGE ALL ? (NO):')
CALL CMREAD
IF(NOCMN.EQ.0)GOTO I27
IF(CMN.EQ.'X') RETURN
IF((CMN.EQ.'Y').OR.(CMN.EQ.'YES'))GOTO 200

127 WRITE(JLUN,128)
128 FORMAT(/,'ENTER FUNCTION INDEX :')
CALL CMREAD

IF(CMN.EQ.'X')RETURN
IF(NOCMN.EQ.0)GOTO 127
DECODE(80.*.CMN.ERRIIAO)MFCN
IF (MFCN.LT.0) THEN
PRINT A,'ERROR] INCORRECT FUNCTION NUMBER'
GOTO I27
ELSE IF(MFCN.EQ.O) THEN
RETURN
ENDIF

DD I30 JJI1,NEQS
lF(|EQ2$(JJ,l) .EQ. MFCN) THEN
NOF-JJ
GOTO I50
ENDIF
130 CONTINUE

1A0 PRINT *,'ERROR] BAD FUNCTION LABEL'
GOTO.127

150
155

I60
165

C

C

WRITE(JLUN,155)lEQ2$(NOF,l)
FORMAT(/,'FUNCTION ',12,' NOW DEFINED AS :')
CALL FCNWTR(NOF)

CALL EQRVS(NOF)
IF(OPT2.EQ.'X')THEN
DPT2-' '

RETURN
ENDIF

WRITE(JLUN,165)

FORMAT(/.'FUNCTION TYPE 2')

CALL CMREAD

IF(NOCMN.EQ.0)THEN
CMN=FCNSET(IFCNTP(NOF))
FPOCMN=I
LPOCMN=6

ENDIF

CALL CMNCPR(FCNSET,30)

IF (IFOUND.EQ.0)THEN
WRITE(JLUN,900I)CMN
GOTO I60

ELSE IF(IFOUND.GT.1)THEN
WRITE(JLUN,9002)(CMATCH(K),K=1,IFOUND)
GOTO 160

ENDIF

CALL FCNLIB(NOF)
IF(OPT2.EQ.'N')THEN
OPT2I'U'
RVS(NOF)=0.
GOTO I50
ELSE IF(OPT2.EQ.'Q')THEN
OPT2I'U'
RETURN
ENDIF

GOTO A00

C-~-CHANGE ALL

200

C
210
215

220
225

CONTINUE
DO 300 NOF-I,NEQS

WRITE(JLUN,215)IEQ2$(NOF,1)
FORMAT(/.'FUNCTION ',12.' NOW DEFINED AS :')
CALL FCNWTR(NOF)

CALL EQRVS(NOF)
IF(OPT2.EQ.'X‘)THEN
DPT2=' '

RETURN
ENDIF

WRITE(JLUN,225)

FORMAT(/,'FUNCTION TYPE ?')

CALL CMREAD

C 2A

IF(NOCMN.EQ.0)CMN=FCNSET(IFCNTP(NOF))
CALL CMNCPR(FCNSET,30)
IF (IFOUND.EQ.0)THEN
WRITE(JLUN,9OOI)CMN
GOTO 220
ELSE IF(IFOUND.GT.I)THEN
WRITE(JLUN,9002)(CMATCH(K),K=1,IFOUND)

GOTO 220
ENDIF
C
CALL FCNLIB(NOF)
IF(OPT2.EQ.'N')THEN
OPTZI' '
RVS(NOF)-0.
GOTO 210
ELSE IF(OPT2.EQ.'Q')THEN
OPT2I'U'
RETURN
ENDIF
C
300 CONTINUE
C
RETURN
c ..........
A00 CONTINUE
C

A05 WRITE(JLUN,AI0)

AIO FORMAT(/,'CHANGE ANOTHER FUNCTION ? (YES):')
CALL CMREAD
IF(NOCMN.EQ.0)GOTO I27
|F((CMN.EQ.'Y').OR.(CMN.EQ.'YES'))GOTO I27

RETURN
C ----- FORMAT -----
C
900I FORMAT(/,' ERROR J NO SUCH FUNCTION : ',A80)

9002 FORMAT(/.' ERROR J FUNCTIONS MIXED UP'./.(A6,' ?'))

END

C

C---EQRVS ----------------------------------------------------
C

SUBROUTINE EQRVS(K)

C

C ----- WRITTEN BY : TSUN-YU CHOU, SEPT. I982

c

$INSERT COMFILE

$INSERT COMAREA

c ............................................................
IF(IPNT(I).EQ.I)WRITE(JLUN,I50)RVS(K)
WRITE(JLUN,101)

lOl FORMAT(/.'REVERSE INPUT AND OUTPUT ? (NO):')

CALL CMREAD

IF(NOCMN.EQ.0)RETURN
IF((CMN.EQ.'Y').OR.(CMN.EQ.'YES'))GOTO 100
IF(CMN.EQ.'X')OPT2I'X'

RETURN

100 WRITE(JLUN,110) IEQZS(K,1)
llO FORMAT(/.'FUNCTION ',IZ.‘ NOW DEFINED AS :')

RR-RVS(K)

OPT2I'U'
IF(RR.EQ.0.)THEN
RVS(K)=I.

ELSE IF(RR.EQ.I)THEN
RVS(KlIO.

ENDIF

CALL FCNWTR(K)
150 FORMAT(/,'*** RVS I ',I1,' ***')

C
RETURN
END
C
C---SHOW ---------------------------------------------------
C
SUBROUTINE SHOW
C
C ----- WRITTEN BY : TSUN-YU CHOU, JULY 1982
C
C DISPLAY USER-DEFINED FUNCTIONS
C AND SYSTEM-FORMATED FUNCTIONS
C

$INSERT COMFILE

$INSERT COMAREA

c ............................................................
OPTII'LOOK '

A25 WRITE(JLUN.A26)
A26 FORMAT(/.'LOOK AT ALL ? (NO):')
CALL CMREAD
IF(NOCMN.EQ.0)GOTOA3O
IF(CMN.EQ.'X') RETURN
IF((CMN.EQ.'Y').OR.(CMN.EQ.'YES'))GOTO 560

A30 WRITE(JLUN,A32)
A32 FORMAT(/.'ENTER FUNCTION INDEX :')
CALL CMREAD

IF(CMN.EQ.'X')RETURN
IF(NOCMN.EQ.0)GOTO A30
DECODE(80.*.CMN,ERRI5A0)MFCN
IF (MFCN.LT.O) THEN
PRINT *,'ERROR] INCORRECT FUNCTION NUMBER'
GOTO A30

ELSE IF(MFCN.EQ.O) THEN
RETURN
ENDIF

525 OD 530 JJI1,NEQS
IF(IEQ25(JJ,1) .EQ. MFCN) THEN
NOFIJJ
GOTO 550
ENDIF
530 CONTINUE
5A0 PRINT *,'ERROR] BAD FUNCTION LABEL'
GOTO A30

550 WRITE(JLUN,9003) IEQ2$(NOF,1)
CALL FCNLOOK(NOF)
GOTO 600

C ----- LOOK AT ALL FUNCTIONS.

560 CONTINUE
DO 570 NOF=I,NEQS
WRITE(JLUN,9003) lEQ2$(NOF,l)
CALL FCNLOOK(NOF)

570 CONTINUE

C

RETURN
c ..........
600 CONTINUE
C

605 WRITE(JLUN,6lO) ,

610 FORMAT(/,'LOOK AT ANOTHER FUNCTION ? (YES):')
CALL CMREAD
IF(NOCMN.EQ.0)GOTO A30
IF((CMN.EQ.'Y').OR.(CMN.EQ.'YES'))GOTO A30

RETURN
C
C ----- FORMAT -----
C
9001 F0RMAT(/.' ERROR ] N0 SUCH FUNCTION : ‘,A80)

9002 FORMAT(/.' ERROR J FUNCTIONS MIXED UP',/,(A6,' ?'))
9003 FORMAT(//,' FUNCTION ',IZ,’ :')
C

END
C

C

C---FCNWTR ----------------------------------------------------
C

SUBROUTINE FCNWTR(K)

c

C---WRITTEN BY : TSUN-YU CHOU. AUGUST I982

C

$INSERT COMFILE
$INSERT COMAREA
c ............................................................

NNEQTPINEQTP(K)
C
IF(NNEQTP .LE. 6) THEN
IF(RVS(K).EQ.0)THEN
IF(OPTI.EQ.'DEFINE')THEN
WRITE(JLUN,9OIO)CHIO(K,1),NFCN(K),NFCN(K),
+ CHIO(K,2),NFCN(K)
ELSE IF(OPTl.EQ.'CHANGE')THEN
WRITE(JLUN,9OZO)CHIO(K,1),NFCN(K),FCNSET(IFCNTP(K)),
+ CHIO(K,2),NFCN(K)
ENDIF
C
ELSE IF(RVS(K).EQ.I)THEN
WRITE(JLUN,9020)CH|O(K.2),NFCN(K),FCNSET(IFCNTP(K)),
+ CHIO(K,1),NFCN(K)
ENDIF
C
c
C
ELSE IF(NNEQTP .GE. I2) THEN
IF(RVS(K).EQ.O) THEN
IF(OPT1.EQ.'DEFINE')THEN
WRITE(JLUN.8810)CHIO(K,1),IEQZS(K,1),IEQZS(K,I),
+ CHIO(K,2),IEQ25(K,2),
+ CHIO(K,3),IEQ2$(K,2),IEQ2$(K,I),
+ CHIO(K,A),IEQ25(K,1)
ELSE IF(OPTl.EQ.'CHANGE')THEN
WRITE(JLUN,882O)CHIO(K,I),IEQ25(K,I),FCNSET(IFCNTP(K)),
+ CHIO(K,2),IEQZS(K,2),
+ CHIO(K,3),IEQ25(K,2),FCNSET(IFCNTP(K)),
+ CHIO(K.A),IEQ2$(K,I)
ENDIF
C
ELSE IF(RVS(K).EQ.I) THEN
WRITE(JLUN,8820)CHIO(K,2),IEQZS(K,2),FCNSET(IFCNTP(K)),
+ CHIO(K,l),IEQ2$(K,l),
+ CHIO(K,A),IEQ2$(K,I),FCNSET(IFCNTP(K)),
+ CHIO(K.3).IEQZS(K,2)
C
ENDIF
c
ENDIF
c
c ----- FORMAT ----------
C

9010 FORMAT(/.5X,A1,12,'
9020 FORMAT(/.5X,A1,12,'
8810 FORMAT(/.5X.AI.I2.'

FCN',I2,' (',AI.I2.')' )
',A6.' ('.AI.I2.')')
FCN',12,' (',Al.|2.')'./.

+ 5X,Al,|2,' FCN',12,' (',A1,|2,')' )
8820 FORMAT(/,5X,A1,I2,' ',A6,‘ (',A1,I2,')',
+ /’5x9A]’|29' ',A69' (',A‘Dlzi')')

C
RETURN

END
C
C---FCNLIB ---------------------------------------------------
C

SUBROUTINE FCNLIB(I)

C
C ----- WRITTEN BY : TSUN-YU CHOU. SEPT. I982
C
$INSERT COMFILE
$INSERT COMAREA
C
C
C VARIABLE : IFCNTP --- STORE FUNCTION TYPE BELOW.
C
C AAAAA ALL FUNCTION INDICES USED IN OTHER SUBROUTINES ARE
C AAAAA ACCORDING TO THE FOLLOWING TYPE NUMBER.
C
C TYPE DEFINITION
C ............................................................
C 01 CONST Y-CO
C 02 LINEAR Y-C0+CIAX
C 03 SIN Y-CIASIN(C2AX+C3) + CA
C 0A COS Y-CIACOS(C2AX+C3) + CA
C 05 ASIN Y-CIAASIN(C2AX+C3) + CA
C
C 06 ACOS Y-CIAACOS(C2AX+C3) + CA
C 07 A852 Y-CIAXAABS(X)
C 08 ABSQRT Y-CIASIGN(X)ASQRT(ABS(X))
C 09 EXP Y-CIAEXP(C2AX+C3) + CA
C IO ALOG Y=CIAALOG(C2AX+C3) + CA
C
C II DSIN Y-CI/SIN(C2AX+C3) + CA
C I2 DCOS Y-CI/COS(C2AX+C3) + CA
C I3 DASIN Y-CI/ASIN(C2AX+C3) + CA
C IA DACOS Y-CI/ACDS(C2AX+C3) + CA
C I5 POLY Y-CO+CIAX+C2AXAA2+C3AXAA3+CAAXAAA
C 16 COLUMB Y=CIA SIGN(X)
C
C 17 BRKPT Y-YI+CIA(X-XI) , X <=x1
C Y-YI+C2A(X-XI) , XI< X
C
C 18 SATUR YIYl , x <- XI
C Y=YI+((Y2-YI)/(x2-XI))A(X-XI) , XI<X<X2
C YIY2 , X2 <= X
C
C I9 ADJSAT Y-YI+CIA(X-XI) . x <- XI
C YIYI+((Y2-Yl)/(X2-Xl))*(X-Xl) . XI<X<X2
C Y-Y2+C2A(X-X2) , X2<- X
C
C 20 PWLIN PIECEWISE LINEAR
C
C 2I INTPLO INTERPOLATION POLYNOMIAL
C

29

Y=Co + CI/(I+(X/C2)AA2)

RCPSQR

22

CV

5
N").
IOT
I_T.|.
WEU
P50.
TAO)
A22
500
JTT
nuF.F.
.A.5.5

RCPSQR
SET027

DATA FCNSET/

+

+++

Pu

DATA ((PARDF(I,J),JII,10),III,30)

momoo
EEEEE
EEEEE
EEEEE
EEEEE
EEEEE
EEEEE
EEEEE
EEEAE
+ + + + +

C--11

111]]

+ + + + +

C--16

omoqd
EEEEE
EEEEE
EEEEE
EEEEE

,,,,

”””

+ +

C--21

”””

+ + + + +

C--26

0
0

0

o

0

A
A
m
m

0

+

991/
Emma
«mum
mwmm
mama
mmmm
mmmm

0000

””
0000
0000

0000

+ + + +

C---_-_--------_------_---_----------------------------------

.EQ.'CONST')THEN

IF(CMATCH(IFOUND)

IFCNTP(I)I1

CALL FCNSUBI(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)IZ
CALL FCNSU82(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)I3
CALL FCNSUB3(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)IA
CALL FCNSUBA(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)I5
CALL FCNSU85(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)I6
CALL FCNSUB6(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)I7
CALL FCNSUB7(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)-8
CALL FCNSU88(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)-9
CALL FCNSU89(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)IIO
CALL FCNSUBI0(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)III
CALL FCNSUBII(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)IIZ
CALL FCNSUB12(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)I13
CALL FCNSUBI3(I)

ELSE IF(CMATCH(IFOUND)
IFCNTP(I)IIA
CALL FCNSUBIA(I)

ELSE IF(CMATCH(IFOUND)

.EQ

.EQ

.EQ

.EQ

.EQ

.EQ

.EQ

.EQ

.EQ

.EQ

.EQ

.EQ

.EQ

.EQ

C 30

.'LINEAR')THEN

.'SIN')THEN

.'COS')THEN

.'ASIN')THEN

.'ACOS')THEN

.'ABSZ')THEN

.'ABSQRT')THEN

.'EXP')THEN

.'ALOG')THEN

.'DSIN')THEN

.'DCOS')THEN

.'DASIN')THEN

.'DACOS')THEN

.'POLY')THEN

IFCNTP(I)I15
CALL FCNSU815(I)

C
ELSE IF(CMATCH(IFOUND) .EQ.'COULOM')THEN
IFCNTP(I)II6
CALL FCNSUBI6(I)
C
ELSE IF(CMATCH(IFOUND) .EQ.'BRKPT')THEN
IFCNTP(I)=I7
CALL FCNSUBI7(I)
C
ELSE IF(CMATCH(IFOUND) .EQ.'SATUR')THEN
IFCNTP(I)IIB
CALL FCNSUBI8(I)
C
ELSE IF(CMATCH(IFOUND) .EQ.'ADJSAT')THEN
IFCNTP(I)II9
CALL FCNSUBI9(I)
C
ELSE IF(CMATCH(IFOUND) .EQ.'PWLIN')THEN
IFCNTP(I)IZO
OPT3I'LINEAR'
CALL FCNSUBZO(I)
C
ELSE IF(CMATCH(IFOUND) .EQ.'|NTPLO')THEN
IFCNTP(I)IZI
OPT3-'INTPLO'
CALL FCNSUBZO(I)
C
ELSE IF(CMATCH(IFOUND) .EQ.'RCPSQR')THEN
IFCNTP(I)I22
CALL FCNSU822(I)
C
ELSE IF(CMATCH(IFOUND) .EQ. 'QUIT')THEN
OPT2-'Q'
RETURN
C
ENDIF
C
RETURN
END
C
C
C---FCNLOOK --------------------------------------------------
C
SUBROUTINE FCNLOOK(I)
C
C ----- WRITTEN BY : TSUN-YU CHOU. SEPT. I982
C

$INSERT COMFILE

$INSERT COMAREA

C ............................................................
C

IF(IFCNTP(I).EQ.I) THEN

CALL FCNSUBI(|)

ELSE IF(IFCNTP(I) EQ.

CALL FCNSUBZ(|)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUB3(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUBA(I)

2)

3)

A)

ELSE IF(IFCNTP(I).EQ.5)

CALL FCNSUB5(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUB6(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUB7(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUBB(|)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUB9(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUBIO(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUB11(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUB12(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUBI3(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUBIA(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUB15(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUBI6(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSUBI7(I)

ELSE IF(IFCNTP(I).EQ.

CALL FCNSU818(I)

6)

7)

8)

9)

IO)

II)

12)

13)

IA)

15)

l6)

17)

I8)

THEN

THEN

THEN

THEN

THEN

THEN

THEN

THEN

THEN

THEN

THEN

THEN

THEN

THEN

THEN

THEN

THEN

32

nonnnn

C
ELSE IF(IFCNTP(I).EQ.19) THEN
CALL FCNSUB19(I)
C
ELSE IF(IFCNTP(I).EQ.20) THEN
CALL FCNSUB20(I)
C
ELSE IF(IFCNTP(I).EQ.2I) THEN
CALL FCNSUBZO(I)
C
ELSE IF(IFCNTP(I).EQ.22) THEN
CALL FCNSU822(I)
C
ENDIF
C
RETURN
END
C
C---PARTAK --------------------------------------------------
C
SUBROUTINE PARTAK(LIBN,K)
C
c ----- WRITTEN BY: TSUN-YU CHOU, OCT. 1982
SUBROUTINE 'PARTAK' AND 'SUBWTR' CALLED BY FCNSUB'S
(3, A, 5, 6, 9, IO, 11, 12, 13. IA)
FOR ENTER PARAMETERS AND VERIFY THE VALUES.
ARGUMENTS : LIBN---ROW NUMBER OF 'PARDF'.(FUNCTION TYPE INDEX)
K ----- EQUATION INDEX. EX: FCN 2<--(2 IS INDEX)
C

$INSERT COMFILE
$INSERT COMAREA
C ............................................................
C
10 NP-I
100 WRITE(JLUN,IOOZ)PARDF(LIBN,1)
1002 FORMAT(/.'VALUE OF CI ? (',F8.A,' ):')
CALL CMREAD
IF(NOCMN.EQ.0)THEN
CI-PARDF(LIBN,I)
GOTO 15
ENDIF
DECODE(80.A.CMN.ERR=IIO)CONSTI
C1=CONSTI
I5 IF(CI.EQ.O) THEN
CALL ERRMSG(s)
GOTO 10
ENDIF

20 NPI2

200 WRITE (JLUN,2002) PARDF (LIBN,2)

2002 FORMAT(‘VALUE OF C2 7 ('.F8.A,' ):')
CALL CMREAD

25

30
300
3002

A0
A00
A002

IF(NOCMN.EQ.0)THEN
C2-PARDF(LIBN,2)
GOTO 25
ENDIF
DECODE(80,A,CMN,ERR=IIO)CONST2
C2-C0NST2
IF(C2 EQ.0) THEN
CALL ERRMSG(5)
GOTO 20
ENDIF

NP-3
WRITE(JLUN,3OOZ)PARDF(LIBN,3)
FORMAT(‘VALUE OF C3 ? (',F8.A,' ):')
CALL CMREAD
IF(NOCMN.EQ.0)THEN

C3-PARDF(LIBN,3)

GOTO A0
ENDIF
DECODE(80,A,CMN,ERR=IIO)CONST3
C3-CONST3

NP=A
WRITE(JLUN.A002)PARDF(LIBN.A)
FORMAT(‘VALUE OF CA ? ('.F8.A.' ):')
CALL CMREAD
IF(NOCMN.EQ.0)THEN

CA=PARDF(LIBN,A)

GOTO 50
ENDIF
DECODE(80.A,CMN,ERR=IIO)CONSTA
CA=CONSTA

CONTINUE
IF(RVS(K).EQ.I) THEN
IF(LIBN.LE.IO) THEN
EQPAR(K,I)-I./C2
EQPAR(K,2)-I./CI
EQPAR(K,3)--CA/CI
EQPAR(K,A)--C3/C2
ELSE IF((LIBN.GE.II).AND.(LIBN.LE.IA))THEN
EQPAR(K.I)-I./C2
EQPAR(K,2)-C1
EQPAR(K.3)ICA
EQPAR(K,A)I-C3/C2
ENDIF
ELSE
EQPAR(K,I)-CI
EQPAR(K,2)-C2
EQPAR(K,3)-C3
EQPAR(K,A)ICA
ENDIF

RETURN

110 CALL ERRMSG(A)
GOTO(10,20,30,A0),NP

C
END
C .
C---SUBWTR --------------------------------------------------
C
SUBROUTINE SUBWTR(LIBN.K)
C

$INSERT COMFILE

$INSERT COMAREA
CHARACTERA6 FTYPE(1A),RFTYPE(IA)
COMMON /SP/ FTYPE,RFTYPE

DATA RFTYPE /' '.' ','A ASIN'.'A ACOS',
+ l* SINI,I* COSI’I I’I I,
+ 'A ALOG','A EXP ','/ ASIN','/ ACOS’,
+ '/ SIN ','/ COS '/
c _
C ............................................................
C
60 IF(RVS(K).EQ.I) THEN
C
IF(LIBN.LE.IO)THEN
TTl=l./EQPAR(K,2)
TT2=l./EQPAR(K,1)
TT3=-EQPAR(K,A)ATT2
TTAI-EQPAR(K,3)*TTI
ELSE IF(LIBN.GE.11)THEN
TTIIEQPAR(K,2)
TT2=I./EQPAR(K,I)
‘ TT3--EQPAR(K,A)ATT2
TTA-EQPAR(K,3)
ENDIF
C
IF(OPT2.EQ.'U') THEN
WRITE(JLUN.1003)CH|O(K,2),IEQZS(K,2),TT1,
+ FTYPE(IFCNTP(K)).TT2,
+ CHIO(K,1),IEQZS(K,1),TT3,TTA
C
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,IOO3)CHIO(K,A),IEQ25(K,I),TT1,
+ FTYPE(IFCNTP(K)).TT2,
+ CHIO(K.3).IEQ2$(K,2),TT3,TTA
ENDIF
C
ELSE IF(OPT2.EQ.'S') THEN
WRITE(JLUN,1005)CHIO(K,1),IEQZS(K,l),EQPAR(K,1),
+ RFTYPE(IFCNTP(K)).EQPAR(K.2),
+ CHIO(K,2),IE025(K,2),EQPAR(K.3).EQPAR(K,A)
C

IF(NOIO(K) EQ.A) THEN

C 36

WRITE(JLUN.1005)CHIO(K,3),IEQ2$(K,2),EQPAR(K,I),

+ RFTYPE(IFCNTP(K)).
+ EQPAR(K,2),CH|O(K,A),IEQ2$(K,1),EQPAR(K,3),
+ EQPAR(K,A)
ENDIF
C
ENDIF
ELSE
WRITE(JLUN,lOO3)CHIO(K,l),IEQZS(K,1),EQPAR(K,1),
+ FTYPE(IFCNTP(K)).EQPAR(K,2).
+ CHIO(K,2),IEQZS(K,2),EQPAR(K.3).EQPAR(K,A)
C
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,lOO3)CHIO(K,3),IEQZS(K,2),EQPAR(K,l),
+ FTYPE(IFCNTP(K)),
+ EQPAR(K,2),CHIO(K,A),IEQ2$(K,I).EQPAR(K,3),
+ EQPAR(K,A)
ENDIF
c
ENDIF
C
1003 FORMAT(/.5X.AI,12,' -',1PEI2.A,' ',A6,'(',lPE12.A,
+ ' A ',Al,|2,' +',1PE12.A,' )',/,9X,'+',1PE12.A)
1005 FORMAT(/,5X,Al,12,' I',1PE12.A,' ',A6,'(',1PE12.A,
+ ' A ',Al,12,' +',lPE12.A,' )',/,9X,'+',1PE12.A)
C
RETURN
END
C
C---IORVSR --------------------------------------------------
C
SUBROUTINE IORVSR(K.15)
C
C ----- WRITTEN BY: TSUN-YU CHOU, DEC. 1982
C ARGUMENTS : K ----EQUATION INDEX.
C 15----FUNCTION TYPE INDEX.
C

$INSERT COMFILE

$INSERT COMAREA

C

CHARACTER*6 FTYPE(IA),RFTYPE(1A)

COMMON /SP/ FTYPE,RFTYPE

DATA FTYPE/' ',' ','* SIN ','* COS ',
1* ASINI’I* ACOSI,I I’l I,
'* EXP ','* ALOG','/ SIN ','/ COS ',

+ '/ ASIN','/ ACOS'/

+
+

c ............................................................
IF(RVS(K).EQ.I)THEN
11:2
12-1
I3=A
IA=3

ELSE
IlII
12I2
|3=3
IAIA

ENDIF

IF(I5.EQ.7)THEN
WRITE(JLUN,IOOZ)CHIO(K,Il),IEQZS(K,II),CHIO(K,|2),
+ IEQ2$(K,12),CHIO(K.12),IEst(K,12)
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,1002)CHIO(K,I3).IEQZS(K,I2),CHIO(K,IA),
+ IEQ2$(K,Il),CHIO(K,IA),IEQZS(K,II)
ENDIF

ELSE IF(15.EQ.8)THEN
WRITE(JLUN,1003)CHIO(K,ll),IEQZS(K,II),CHIO(K,|2),
+ IEQZS(K,12),CH|O(K,12),IEQZS(K,|2)
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,lOO3)CH|O(K,I3),IEQZS(K,12),CHIO(K.IA).
+ IEQ2$(K,Il),CHIO(K,IA),IEQ2$(K,II)
ENDIF

ELSE
WRITE(JLUN,1001)CHIO(K,Il),IEQZS(K,Il),FTYPE(15),

+ CHIO(K,I2).IEQ2$(K,I2)
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,IOOl)CHIO(K,I3).IEQZS(K,12),FTYPE(15),

+ CHIO(K,IA),IEQ2$(K,II)

ENDIF
C
ENDIF
C
1001 FORMAT(/,5X,A1,12,' = C1 ',A6,' ( C2 A ',Al,12,

+ ' + C3 ) + CA')

1002 FORMAT(/.5X,A1,IZ,' . C1 A ',Al,12,' A ABS( ',Al,|2,' )')
1003 FORMAT(/,5X,Al,|2,' = C1 A SIGN( ',Al,12,' ) A SQRT( ',

+ 'ABS( ',Al.12.' ) )' )

C

RETURN

END
C
C---GENPAR --------------------------------------------------
C

SUBROUTINE GENPAR(K,KDF.KPD.NCOOR)
C THIS SUBROUTINE USED FOR ENTER PAIR OF x.Y COORDINATE.
C ARGUMENTS : K ----- EQUATION INDEX.
C KDF---- ROW NUMBER OF 'PARDF'.
C KP0--—- COLUMN NUMBER OF 'PARDF'.
C NCOOR-- COORDINATE INDEX. EX: (XI,Y1)

$INSERT COMFILE
$INSERT COMAREA
C ............................................................

C 38

IO NPIl
WRITE(JLUN,1002)NCOOR,PARDF(KDF,KPO)
1002 FORMAT(/,'VALUE OF X',Il,' ? (',F8.A,' ):')
CALL CMREAD
IF(NOCMN.EQ.0)THEN
EQPAR(K.KPO)-PARDF(KDF.KP0)
GOTO 20
ENDIF
DECODE(80,A,CMN,ERR=IIO)CONSTI
EQPAR(K,KPO)ICONST1

20 NPI2
KPoz-KPO+I
WRITE(JLUN,2002)NCOOR.PARDF(KDF,KP02)
2002 FORMAT(‘VALUE OF Y',|l,' 2 (',F8.A,' ):')
CALL CMREAD
IF(NOCMN.EQ.0)THEN
EQPAR(K,KP02)IPARDF(KDF,KP02)
GOTO 30
ENDIF
DECODE(80.A,CMN,ERR-IIO)CONST2
EQPAR(K,KPOZ)ICONST2

C
30 CONTINUE
C
RETURN
C

110 CALL ERRMSG(A)
GOTO (IO.20).NP

C
END

C

C---GENCON --------------------------------------------------

C
SUBROUTINE GENCON(K.KDF.KPO.NCONS)

C THIS SUBROUTINE USED FOR ENTER CONSTANT Ci.

C ARGUMENTS : K ----- EQUATION INDEX.

C KDF---- ROW NUMBER OF 'PARDF'.

C KPO---— COLUMN NUMBER OF 'PARDF'.

C NCONS—--CONSTANT INDEX. Ex: C1, C2

$INSERT COMFILE
$INSERT COMAREA
C ............................................................
10 WRITE(JLUN,1002)NCONS,PARDF(KDF,KPO)
1002 FORMAT(/,'VALUE OF C',|l,' ? (',F8.A,' ):')
CALL CMREAD
IF(NOCMN.EQ.0)THEN
EQPAR(K,KPO)-PARDF(KDF,KPO)
GOTO 20
ENDIF
DECODE(80.A.CMN,ERR=110)CONSTI
EQPAR(K,KPO)ICONST1

20 CONTINUE

RETURN
C
110 CALL ERRMSG(A)
GOTO 10
C
END
C
C---COOREX --------------------------------------------------
C
SUBROUTINE COOREX(K.IPO,JPO)
C THIS SUBROUTINE EXCHANGE THE X—Y COORDINATE.
C ARGUMENTS : K ---:- EQUATION INDEX.
C IPO--—- STORAGE POSITION OF x COORDINATE.
C JPO---- STORAGE POSITION OF Y COORDINATE.
$INSERT COMAREA
c ............................................................

TEMPI=EQPAR(K.IPO)
EQPAR(K,IPO)IEQPAR(K,JPO)
EQPAR(K.JP0)ITEMPI

C
RETURN
END
C
C---FCNSUBI --------------------------------------------------
C
SUBROUTINE FCNSUBI(K)
C
C ----- FUNCTION TYPE : CONST Y=CO.
C
C ----- WRITTEN BY : TSUN-YU CHOU, AUGUST 1982
C ----- MODIFIED : FEB. 1983.
C

$INSERT COMFILE
$INSERT COMAREA

C ............................................................
IF(OPTl.EQ.'LOOK ') GDTO 200
C
IF(RVS(K).EQ.I)THEN
OPT2-'N'
CALL ERRMSG(6)
RETURN
ENDIF
C
WRITE(JLUN,lOOl)CHIO(K,l),IEQ25(K,1)
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,lOOl)CHIO(K,3),IEQ25(K,2)
ENDIF
C
1001 FORMAT(/.5X,Al,12,' = C0 ')
C

100 CALL GENCON(K,1,1,0)

IF(EQPAR(K.I).EQ.O)THEN
CALL ERRMSG(5)
GOTO 100

ENDIF

CALL ERRHSG(7)

CALL CMREAD

IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 200
RETURN

200 WRITE(JLUN.IOO3)CHIO(K,1),IEQ2$(K,I),EQPAR(K,I)
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,lOO3)CHIO(K,3),IEQ2$(K,2),EQPAR(K,1)

ENDIF

C
I003 FORMAT(/.5X,Al,l2,' =',lPE12.A)
C

RETURN

END
C
C---FCNSUBZ --------------------------------------------------
C

SUBROUTINE FCNSU32(K)
C
C ----- FUNCTION TYPE : LINEAR Y=C0 +CIAX

$INSERT COMFILE
$INSERT COMAREA

c ............................................................
IF(OPT1.EQ.'LOOK ') GOTO 60
C
IF(RVS(K).EQ.I)THEN
11-2
12-1
I3-A
lA-3
ELSE
11-1
12-2
)3-3
IA-A
ENDIF
C .
WRITE(JLUN,1001)CHIO(K.Il),IEQ2$(K,Il),CHIO(K,12),IEQZS(K,12)
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,1001)CHIO(K,I3).IEQZS(K,12),CHIO(K,IA),IEQZS(K,Il)
ENDIF
C
1001 FORMAT(/.5X,Al,12,' = C0 + C1 A ',Al,12)
C
CALL GENCON(K,2,1,0)
C

10 CALL GENCON(K.2,2,I)
15 IF(EQPAR(K,2).EQ.O) THEN

CALL ERRMSG (5)

GOTO 10
ENDIF
C
IF(RVS(K).EQ.I) THEN
EQPAR(K,2)=I./EQPAR(K,2)
EQPAR(K,1)=-EQPAR(K,I)AEQPAR(K,2)
ENDIF
C
CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 60
RETURN
C
60 1F(RVS(K).EQ.I .AND. OPT2.EQ.'U')THEN
TT2-1./EQPAR(K.2)
TT1--EQPAR(K,I)ATT2
C
WRITE(JLUN.IOO3)CHIO(K,2),IEst(K,2),TTI,TT2.
+ CH10(K,I),IEQ2$(K,I)
C
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,IOO3)CHIO(K,A),IEQ2$(K,1),TT1,TT2,
+ CHIO(K.3).IEQ2$(K,2)
ENDIF
C
ELSE
WRITE(JLUN,1003)CHIO(K,1),IEQZS(K,1),EQPAR(K,I),
+ EQPAR(K,2),CHIO(K,2),IEQZS(K,2)
C
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN.lOO3)CHIO(K,3),IEQZS(K,2),EQPAR(K,1),
+ EQPAR(K,2),CHIO(K.A),IEQ2$(K,1)
ENDIF
C
ENDIF
C
1003 FORMAT(/,5X,AI,12.' =',lPE12.A,' + ',IPE12.A,' A ',Al,12)
RETURN
END
C
C---FCNSUB3 --------------------------------------------------
C
SUBROUTINE FCNSUB3(K)
C
C ----- FUNCTION TYPE : SIN
C Y-CIASIN(C2AX+C3) + CA

$INSERT COMFILE

$INSERT COMAREA

C ............................................................
IF(OPT1.EQ.'LOOK ') GOTO 60
CALL IORVSR(K.3)
CALL PARTAK(3.K)

CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 6O

RETURN
C
60 CALL SUBWTR(3,K)
RETURN
END
C
C---FCNSUBA --------------------------------------------------
C
SUBROUTINE FCNSUBA(K)
C
C ----- FUNCTION TYPE : COS
C Y=CIACOS(C2AX +C3) + CA

$INSERT COMFILE
$INSERT COMAREA

c ............................................................
IF(OPTl.EQ.'LOOK ') GOTO 60
CALL IORVSR(K.A)
CALL PARTAK(A,K)
CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 60
RETURN
C
60 CALL SUBWTR(A.K)
RETURN
END
C
C---FCNSUBS --------------------------------------------------
c _
SUBROUTINE FCNSUBS(K)
C
C ----- FUNCTION TYPE : ASIN (SIN INVERSE)
C Y-CIAASIN(C2AX +C3) + CA

$INSERT COMFILE
$INSERT COMAREA

c ............................................................
IF(OPT1.EQ.'LOOK ') GOTO 60
CALL IORVSR(K,5)
CALL PARTAK(5,K)
CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 60
RETURN
C
60 CALL SUBWTR(5.K)
RETURN
END
C
C---FCNSUB6 --------------------------------------------------

C

SUBROUTINE FCNSUB6(K)

C
C ----- FUNCTION TYPE : ACOS
C YICI*ACOS(C2*X +C3) +CA

$INSERT COMFILE
$INSERT COMAREA

C ............................................................
IF(OPTl.EQ.’LOOK ') GOTO 60
CALL IORVSR(K,6)
CALL PARTAK(6,K)
CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 6O
RETURN
C
60 CALL SUBWTR(6,K)
RETURN
END
C
C---FCNSUB7 --------------------------------------------------
C
SUBROUTINE FCNSUB7(K)
C
C ----- FUNCTION TYPE : ABSZ (ABS SQUARE)
C Y=CIAXAABS(X)
C ----- WRITTEN BY : TSUN-YU CHOU. AUGUST 1982
C

$INSERT COMFILE
$INSERT COMAREA
C ............................................................
IF(OPTl.EQ.'LOOK ') GOTO 200
CALL IORVSR(K,7)
IOO CALL GENCON(K,7,1,1)
120 IF(EQPAR(K,I).LE.O)THEN
CALL ERRMSG(5)
GOTO IOO
ENDIF

1F(RVS(K).EQ.1.) EQPAR(K.I)Il./SQRT(EQPAR(K,1))
CALL ERRMSG(7)

CALL CMREAD

IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 200
RETURN

200 IF( (RVS(K).EQ.I).AND.(OPT2.EQ.'U') ) THEN
TT-1./( EQPAR(K,1)AEQPAR(K.I) )

WRITE(JLUN,1003)CHIO(K,2),IEQ2$(K,2),TT,CHIO(K,1),
+ IEQZS(K,1),CHIO(K,1),IEQ25(K,1)
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,1003)CHIO(K,A),IEQZS(K,1),TT,CHIO(K,3),
+ IEQZS(K.2).CHIO(K,3),IEQZS(K,2)
ENDIF

C
ELSE
C
WRITE(JLUN.1005)CHIO(K,1),IEQZS(K,1),EQPAR(K,l),CHIO(K,2),
+ IEQZS(K,2),CHIO(K,2),IEQ2$(K,2)
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,1005)CHIO(K,3),IEQZS(K,1),EQPAR(K,1),CHIO(K,A),
+ IEQZS(K,1),CHIO(K.A),IEQ2S(K,1)
ENDIF
ENDIF
C
1003 FORMAT(/,5X,Al,12,' =‘,lPE12.A,' A ',Al,12,' A ABS( ',
+ ‘ AI.I2.' )' )
1005 FORMAT(/,5X,Al,|2,' =',1PE12.A,' A SIGN( ',Al,12,' )',
+ ' A SQRT( ABS( ',Al,l2,' ) )' )
C
RETURN
END
C
C---FCNSUBB --------------------------------------------------
C
SUBROUTINE FCNSUBB(K)
C
C ----- FUNCTION TYPE : ABSQRT (ABS SQRUARE ROOT)
C Y=C1A SIGN(X) A SQRT( ABS(X) )

$INSERT COMFILE
$INSERT COMAREA
c ............................................................
IF(OPT1.EQ.'LOOK ') GOTO 200
CALL IORVSR(K.8)
100 CALL GENCDN(K,8,I.1)
120 IF(EQPAR(K,I).LE.O)THEN
CALL ERRMSG(5)
GOTO IOO
ENDIF

IF(RVS(K).EQ.I.) EQPAR(K.I)=1./(EQPAR(K.I)AEQPAR(K,1))
CALL ERRMSG(7)

CALL CMREAD

IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 200

RETURN

200 IF( (RVS(K).EQ.I).AND.(OPT2.EQ.'U') ) THEN
TT=1./SQRT(EQPAR(K,1))

C
WRITE(JLUN,lOO3)CHIO(K,2),IEQZS(K,2),TT,CHIO(K,1),
+ IEQ2$(K,I),CHIO(K.I).1Est(K,I)
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,lOO3)CHIO(K,A),IEQZS(K,1),TT,CHIO(K,3),
+ ' IEQ25(K,2),CHIO(K,3),IEQZS(K,2)
ENDIF
C

ELSE

CA5

C
WRITE(JLUN,1005)CHIO(K,I),IEQZS(K,1),EQPAR(K,1),CHIO(K,2).
+ IEQZS(K,2),CHIO(K,2),IEQZS(K,2)
IF(NOIO(K).EQ.A) THEN
WRITE(JLUN,1005)CHIO(K.3).IEQZS(K,1),EQPAR(K,1),CHIO(K,A),
+ IEQZS(K,1),CHIO(K,A),IEQ2$(K,I)
ENDIF
C
ENDIF
C
1003 FORMAT(/,5X,AI,12,"-',IPE12.A,' A SIGN( ',Al,12,' )',
+ 'A SQRT( ABS( ',A1,12,' ) )' )
1005 FORMAT(/,5X,Al,|2,' I',1PE12.A,' A SIGN( ',Al,12,' )',
+ ' A ( '.A1.I2.' ) AA 2')
C
RETURN
END
C
C---FCNSUB9 --------------------------------------------------
C
SUBROUTINE FCNSUB9(K)
C
C ----- FUNCTION TYPE : EXP
C Y=CIAEXP(C2X+C3) + CA

$INSERT COMFILE
$INSERT COMAREA

C ............................................................
IF(OPTI.EQ.'LOOK ') GOTO 60
CALL IORVSR(K,9)
CALL PARTAK(9,K)
CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 60
RETURN
C
60 CALL SUBWTR(9,K)
RETURN
END
C
C---FCNSUBIO --------------------------------------------------
C
SUBROUTINE FCNSUBIO(K)
C
C ----- FUNCTION TYPE : ALOG
C Y-CIAALOG(C2AX+C3) + CA

$INSERT COMFILE
$INSERT COMAREA
C ............................................................
IF(OPT1.EQ.'LOOK ') GOTO 60
CALL IORVSR(K,lO)
CALL PARTAK(IO,K)
CALL ERRMSG(7)
CALL CMREAD

C A6

IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 60

RETURN
C
60 CALL SUBWTR(10,K)
RETURN
END
C
C---FCNSUBII --------------------------------------------------
C
SUBROUTINE FCNSUBII(K)
C
C ----- FUNCTION TYPE : DSIN
C YICI/SIN(C2*X + C3) + CA

$INSERT COMFILE
$INSERT COMAREA

C ............................................................
IF(OPT1.EQ.'LOOK ') GOTO 60
CALL IORVSR(K.11)
CALL PARTAK(11,K)
CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 6O
RETURN
C
60 CALL SUBWTR(II,K)
RETURN
END
C
C--—FCNSUBIz --------------------------------------------------
C
SUBROUTINE FCNSUBIZ(K)
C
C ----- FUNCTION TYPE : DCOS
C Y=CI/COS(C2AX + C3) + CA

$INSERT COMFILE
$INSERT COMAREA

C ............................................................
IF(OPT1.EQ.'LOOK ') GOTO 60
CALL IORVSR(K.12)
CALL PARTAK(12,K)
CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 6O
RETURN
C
60 CALL SUBWTR(12.K)
RETURN
END
C
C---FCNSUBI3 --------------------------------------------------
C

SUBROUTINE FCNSUBI3(K)
C

C ----- FUNCTION TYPE : DASIN

C Y=CI/ASIN(C2AX + C3) + CA
$INSERT COMFILE

$INSERT COMAREA

c ............................................................
IF(OPTI.EQ.'LOOK ') GOTO 60
CALL IORVSR(K.13)
CALL PARTAK(I3,K)
CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 6O
RETURN
C
60 CALL SUBWTR(I3,K)
RETURN
END
C
C---FCNSUB1A --------------------------------------------------
C
SUBROUTINE FCNSUBIA(K)
C
C ----- FUNCTION TYPE : DACOS
C Y-CI/ACOS(C2AX + C3) + CA

$INSERT COMFILE
$INSERT COMAREA

c ............................................................
IF(OPTI.EQ.'LO0K ') GOTO 60
CALL IORVSR(K.1A)
CALL PARTAK(IA,K)
CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOT0 60
RETURN
C
60 CALL SUBWTR(IA,K)
RETURN
END
C
C'I‘FCNSUBIS IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
C
SUBROUTINE FCNSUBI5(K)
C
C """ FUNCTION TYPE : POLY
C Y=C0 + CIAX + CZAXAAZ + C3AXAA3 + CLAXAAL

$INSERT COMFILE
$INSERT COMAREA
C ............................................................
IF(OPT1.EQ.'LOOK ') GOTO 60
IF(RVS(K).EQ.I)THEN
OPT2-'N'
CALL ERRMSG(6)
RETURN
ENDIF

WRITE(JLUN,lOOl)CHIO(K,l),IEQZS(K,1),CH|O(K,2),IEQZS(K,2),
+ CHIO(K,2),IEQZS(K,2).CHIO(K,2),IEQ2$(K,2),
+ CHIO(K.2).IEQ2$(K,2)
lOOl FORMAT(/,5X,Al,12,' - C0 + CIA',AI,12,' + C2A(',AI,12,
+ ')AA2 + C3A('AI,12,')AA3 + CAA(',AI,12.')AAA',/)
C
CALL GENCON(K,15,I,O)
CALL GENCON(K,15,2.I)
CALL GENCON(K,15.3.2)
CALL GENCON(K,15,A,3)
CALL GENCON(K,15,5,A)

CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 60
RETURN
C
60 WRITE(JLUN,1003)CHIO(K,1),IEQ2$(K,I),EQPAR(K,I),EQPAR(K,2),
+ CHIO(K,2),IEQZS(K,2),EQPAR(K,3),
+ CHIO(K,2),IEQZS(K,2),EQPAR(K,A),
+ CHIO(K,2),IEQ2$(K,2),EQPAR(K,5).CH|O(K,2),IEQ25(K,2)
1003 FORMAT(/,5X,Al,12,' -',lPE12.A,' + ',lPE12.A,' A',A1,12,

+ ' +',1PE12.A,' *(',A1,12,')**2',//,9X,
+ '+',IPE12.A,' A(',A1,12,')**3 +',
+ 1PE12.A,' *(',A1,12,')**A')
C
RETURN
END
C
C---FCNSUBI6 --------------------------------------------------
C
SUBROUTINE FCNSUBI6(K)
C
C ----- FUNCTION TYPE : COLUMB
C
C Y I C1 * SIGN(X)
C

$INSERT COMFILE
$INSERT COMAREA

C ............................................................
IF(OPT1.EQ.'LOOK ') GOTO 50
IF(RVS(K).EQ.I)THEN
OPT2='N'
CALL ERRMSG(6)
RETURN
ENDIF
C

WRITE(JLUN,1001)CHIO(K,1),IEQZS(K,1),CHIO(K,2),IEQZS(K,2)
1001 FORMAT(/,5X,AI,12,' I C1 * SIGN(',A1,12,')',/)

. CALL GENCON(K,16,1,I)
CALL ERRMSG(7)

C A9

CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 50
RETURN
C
50 WRITE(JLUN,1003)CHIO(K,1),IEQ2$(K,1),EQPAR(K,1),
+ CHIO(K,2),IEQZS(K,2)
1003 FORMAT(/,5X,Al,|2,' I',1PE12.A,' A SIGN(',A1,|2,')' )

RETURN
END
C---FCNSUBI7 --------------------------------------------------
C
SUBROUTINE FCNSUBI7(K)
C
C ----- FUNCTION TYPE : BRKPT (BREAK POINT)
C Y- YI+C1A(X-x1), WHEN X <=x1
C Y= YI+C2A(x-XI). WHEN x1<x
C
C EQPAR(K,1) STORE XI
C EQPAR(K.2) STORE Y1
C EQPAR(K,3) STORE CI
C EQPAR(K,A) STORE C2
C

$INSERT COMFILE
$INSERT COMAREA
C ............................................................
IF(OPT1.EQ.'LOOK ') GOTO 30
C
WRITE(JLUN,lOOl)CH|O(K,2),IEQZS(K,2),CHIO(K,I),IEQ2$(K,1)
1001 FORMAT(/.5X,'Y= YI+CIA(x-XI), WHEN X <=X1'./.

+ 5x,uy= y1+c2A(X-x1), WHEN XI<X',//.
+ 25X,‘ WHERE X I ',A1,I2./.
'I” 25X,' Y 3',AI,12./)

CALL GENPAR(K,17,I,1)
10 CALL GENCON(K,17,3.1)
IF(RVS(K).EQ.I .AND. EQPAR(K,3).EQ.O)THEN
CALL ERRMSG(B)
GO TO 10
ENDIF
20 CALL GENCON(K,17,A,2)
IF(RVS(K).EQ.I .AND. EQPAR(K,A).EQ.O)THEN
CALL ERRMSG(B)
GO TO 20
ENDIF

IF(RVS(K).EQ.I)THEN
CALL COOREX(K,I,2)
EQPAR(K.3)=I./EQPAR(K,3)
EQPAR(K,A)-I./EQPAR(K,A)
ENDIF

CALL ERRMSG(7)
CALL CMREAD

C 50

IF(CMN.EQ.'Y' .OR. CMN.EQ.‘YES')GOTO 30
RETURN

3O CONTINUE
1F(RVS(K).EQ.I .AND. OPT2.EQ.'U')THEN
TTlIEQPAR(K,2)
TT2=EQPAR(K.I)
TT3=I./EQPAR(K,3)
TTA-I./EQPAR(K,A)
WRITE(JLUN,lOO3)CHIO(K,l),IEQZS(K,1),TT2,TT3,

+ CHIO(K,2),IEQ2$(K,2),TT1,CHIO(K,2),IEQ2$(K,2),
+ TTI,
+ CHIO(K,I),IEQ2$(K,1),TT2,TTA,
+ CHIO(K,2),IEQZS(K,2),TTl,
+ TTI,CH10(K,2),IEQ2$(K,2)
C
ELSE
WRITE(JLUN,IOO3)CHIO(K,I),IEQZS(K,1),EQPAR(K,2),EQPAR(K,3),
+ CHIO(K,2),IEQ2$(K,2),EQPAR(K,1),CHIO(K,2),IEQ2$(K,2),
+ EQPAR(K,I),
+ CHIO(K,I).IEQ25(K,I),EQPAR(K.2),EQPAR(K,A),
+ CHIO(K,2),IEQZS(K,2),EQPAR(K,1),
+ EQPAR(K,I),CHIO(K,2),IEQZS(K,2)
C
ENDIF
C
1003 FORMAT(/,5X,Al,12,' I',1PE12.A,' + ',lPE12.A,' A (',Al,12,
+ ' - ',lPE12.A,' ), ',/,35X,'WHEN ',Al,12,' <- ',lPE12.A,
+ //,5X,Al,|2,' -',IPE12.A,' + ',lPE12.A,' A (',Al,12,
+ ' - ',lPE12.A,' ), ',/.35x,'WHEN ‘,lPE12.A,' < ',Al,12)
C
RETURN
END
C
C---FCNSUBIB --------------------------------------------------
C
SUBROUTINE FCNSUBIB(K)
C
C ----- FUNCTION TYPE : SATUR ( SATURATION )
C
C Y - Yl , X <= x1
C Y - YI+ ((Y2-YI)/(x2-XI))A(X-XI) , x1 < X < X2
C Y - Y2 , X2 <2 X
C
C STORAGE PARDF EQPAR
C 1 x1 X1
C 2 Y1 Y1
C 3 X2 x2
C A Y2 Y2
C

$INSERT COMFILE
$INSERT COMAREA
C ............................................................

IF(OPTl.EQ.'LOOK ') GOTO 60
IF(RVS(K).EQ.I)THEN

OPT2I'N'
CALL ERRMSG(6)
RETURN
ENDIF
C
WRITE(JLUN,lOOl)CHIO(K,2),IEQZS(K,2),CHIO(K,I),IEQZS(K,I)
C
1001 FORMAT(/,5X,'Y - Y1 , x <= Xl',
+ /,5X,'Y - YI+ ((YI-Y2)/(XI—X2))A(x-x1) , X1 < x < X2',
+ /,5X,'Y - Y2 , x2 <= X'//,
+ AOX,’ WHERE X - ',AI,12,/,
+ AOX,‘ Y - ',AI,12,/)
C
CALL GENPAR(K,18,l,l)
CALL GENPAR(K,18,3.2)
C
CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 60
RETURN
C

60 CONTINUE
SLOPE-(EQPAR(K,A)-EQPAR(K,2))/(EQPAR(K,3)—EQPAR(K,I))
WRITE(JLUN,1003)

+ CHIO(K.I),IEQ2$(K,1),EQPAR(K,2).
+ CHIO(K,2),IEQZS(K,2),EQPAR(K,1),
+ CHIO(K,1),IEQ2$(K,1),EQPAR(K,2),SLOPE,
+ CHIO(K,2),IEQZS(K,2),EQPAR(K,1),
+ EQPAR(K,I),CHIO(K,2),IEQZS(K,2),EQPAR(K,3),
+ CHIO(K,1).IEQ2$(K,I),EQPAR(K,A),
+ EQPAR(K,3),CHIO(K,2),IEQZS(K,2)
C
1003 FORMAT(
+ /,5X,Al,12,' I',1PE12.A,' , WHEN ',AI,12,' <I',
+ 1PE12.A,//,
+ 5X,Al,12,' I',1PE12.A,' + ',lPE12.A,' A ( ',AI,12,
+ ' - ',lPE12.A,' )',//,29X,', WHEN ',lPE12.A,‘ < ',AI,12,
+ ' < ',IPE12.A,//,
+ 5X,Al,12,' I',1PE12.A,‘ , WHEN ',lPE12.A,' <=',
+ Al,12)
C
RETURN
END
C
C---FCNSUBl9 --------------------------------------------------
C
SUBROUTINE FCNSUBI9(K)
C
C ----- FUNCTION TYPE : ADJUST
C

C YI Y1 + C1 * (X-XI) , WHEN X <= X1

Y=YI+((Y2-YI)/(x2-XI))A(X-XI) ,

Y: Y2 + C2 A (x-X2)

EQPAR(K.I)
EQPAR(K.2)

EQPAR(K,A)
EQPAR(K,5)
EQPAR(K,6)

STORE
STORE
STORE
STORE
STORE
STORE

VALUE
VALUE
VALUE
VALUE
VALUE
VALUE

OF
OF
OF
OF
OF
OF

X1
Y1
X2
Y2
C1
C2

C 52

WHEN X1 < X < X2
, WHEN X2 <I X

C
C
C
C
C
C EQPAR(K,3)
C
C
C
C

$INSERT COMFILE
$INSERT COMAREA
c ------------------- ,—-----—-—-—-—g --------------------------
IF(OPTl.EQ.'LOOK ') GOTO 80
C
WRITE(JLUN,lOOl)CHIO(K,2),IEQZS(K,2),CHIO(K,1),IEQZS(K,I)
FORMAT(/.5x.'Y= Yl + Cl A (X-Xl) , WHEN X <= Xl',
+ /,5x,'Y=YI+((Y2—YI)/(X2-XI))A(X-x1) , WHEN XI < X < xz',
+ /,5X,'Y= Y2 + C2 A (X-x2) , WHEN X2 <= X ',//,
+
+

1001

AOX,‘ WHERE x = ',Al,IZ,/,
AOX,‘ Y = '.A).I2./)

5 CALL GENPAR(K,19,I,I)
CALL GENPAR(K,19.3.2)
XD=EQPAR(K,3)-EQPAR(K,I)
YDIEQPAR(K,A)-EQPAR(K,2)
1F(RVS(K).EQ.1 .AND. YD.EQ.0)THEN
CALL ERRMSG(B)
GO TO 5
ENDIF
IF(RVS(K).EQ.0 .AND. XD.EQ.0)THEN
CALL ERRMSG(9)
GO TO 5
ENDIF
CALL GENCON(K,19,5,I)
IF(RVS(K).EQ.1 .AND. EQPAR(K,5).EQ.O)THEN
CALL ERRMSG(B)
GO T0 10
ENDIF
CALL GENCON(K.19,6,2)
IF(RVS(K).EQ.I .AND. EQPAR(K,6) EQ.O)THEN
CALL ERRMSG(B)
GO TO 20
ENDIF

10

20

IF(RVS(K).EQ.I)THEN
CALL CDOREX(K.1.2)
CALL COOREx(K,3,A)
EQPAR(K,5)=I./EQPAR(K,5)
EQPAR(K,6)=I./EQPAR(K.6)
ENDIF

CALL ERRMSG(7)

1003

C

+++++++++

+++++++++

+++++++

CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 80
RETURN

CONTINUE
IF(RVS(K).EQ.0)THEN
SLOPE=(EQPAR(K,A)-EQPAR(K,2))/(EQPAR(K,3)-EQPAR(K,1))
ELSE IF(RVS(K).EQ.I)THEN
TTlIEQPAR(K,2)
TT2=EQPAR(K,1)
TT3=EQPAR(K,A)
TTA=EQPAR(K,3)
TT5=I./EQPAR(K,5)
TT6=1./EQPAR(K,6)
SLOPE=(EQPAR(K,3)-EQPAR(K.I))/(EQPAR(K,A)-EQPAR(K,2))

IF(RVS(K).EQ.1 .AND. OPT2.EQ.'U')THEN
WRITE(JLUN,1003)

CHIO(K,2),IEQZS(K,2),TT2,TT5,
CHIO(K,l),IEQ25(K,I),TT1,
CH10(K.1),IEQ2$(K,1),TTI,
CHIO(K,2),IEQZS(K,2),TT2.SLOPE.
CHIO(K.I),1EQ25(K,I),TTI.
TTI,CHIO(K.I),IEQ2$(K.I),TT3,
CHIO(K,2),IEQZS(K,2),TTA,TT6,
CHIO(K.1),|EQ2$(K,I),TT3,
TT3,CHID(K.I),IEQ25(K.I)

WRITE(JLUN,1003)

CHIO(K,l),IEQ2$(K,l),EQPAR(K,2),EQPAR(K,5),
CHIO(K,2),IEQZS(K,2),EQPAR(K,1),
CHIO(K,2),IEQZS(K,2),EQPAR(K,1),
CHIO(K,l),IEQZS(K,l),EQPAR(K.2).SLOPE.
CHIO(K,2),IEQ25(K,2),EQPAR(K,I),
EQPAR(K,1),CHIO(K,2),IEQZS(K,2),EQPAR(K,3),
CHIO(K,1),IEQZS(K,1),EQPAR(K,A),EQPAR(K,6).
CHIO(K,2),IEQ2$(K,2),EQPAR(K,3),
EQPAR(K,3).CHIO(K,2),IEQZS(K,2)

FORMAT(

//,5X,A1,12,' I',IPE12.A,' +',1PE12.A,' A ( ',AI,12,
' - ',1PE12.A,')',//,33X,', WHEN ',AI,12,' <I',1PE12.A,
//,5X,AI,12,' I',1PE12.A,' + ',1PE12.A,' * ( ',AI,12,
' - ',1PE12.A,')',//,33X,', WHEN ',1PE12.A,' < ',AI,12,
' < ',1PE12.A,//,
5X,A1,12,' I',IPE12.A,' +',IPE12.A,' A ( ',AI,12,
' - ',1PE12.A,')',//,33X,', WHEN ',1PE12.A,' <I',AI,I2)

C---FCNSUB20 --------------------------------------------------

C
SUBROUTINE FCNSUBZO(K)
C
C ----- WRITTEN BY: TSUN-YU CHOU, NOV.1982
C ----- FUNCTION TYPE : PWLIN AND INTPLO
C

$INSERT COMFILE
$INSERT COMAREA
c ............................................................
IF(OPTl.EQ.'LOOK ') GOTO 50
IF(RVS(K).EQ.I)THEN
OPT2I'X'
CALL ERRMSG(6)
RETURN
ENDIF

IF(OPT3.EQ.'LINEAR') THEN
ILII
WRITE(JLUN,900)
WRITE(JLUN,920)

ELSE IF(OPT3.EQ.'INTPLO') THEN
ILI2
WRITE(JLUN,910)
WRITE(JLUN,920)

ENDIF

CALL INTERP(K,IL)

CALL ERRMSG(7)

CALL CMREAD

IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOT0 50
RETURN

50 CONTINUE
N-EQPAR(K,1)
IF(IFCNTP(K).EQ.2O) THEN
ILIl
WRITE(JLUN,900)
ELSE 1F(IFCNTP(K).EQ.2I) THEN
1L-2
WRITE(JLUN,910)
ENDIF
WRITE(JLUN,930)N
DO 60 IIl,N
6O WRITE(JLUN,935)I,XCO(IL,I),YCO(|L,I)
IF(IL.EQ.I)THEN
WRITE(JLUN,9AO)EQPAR(K,3),EQPAR(K,A)
ENDIF

C ----- FORMAT -----
900 FORMAT(/.5x.'PIECEWISE LINEAR :' )
910 FORMAT(/.5x,'INTERPOLATION POLYNOMIAL :')
920 FORMAT(/.7X,'MAXIMUM NUMBER OF DATA POINTS IS 20',
+ /.7x.'x COORDINATE MUST BE IN ASCENDING ORDER.')

C 55

930 FORMAT(/.5X,12,' DATA POINTS BEEN INPUT, THERE ARE :',
+ //.I3X,'NO.',10X,'X',I6X,'Y')

935 FORMAT(/,(10X,15,5X,E12.A,5X,E12.A))

9A0 FORMAT(/,6X,'C1I ',F12.A,/,6X,‘C2I ',F12.A)

RETURN
END
C
C
C ............................................................
C
SUBROUTINE INTERP(K,IL)
C
C ----- WRITTEN BY: TSUN-YU CHOU, OCT.I982
C
C NEWTON'S DIVIDED-DIFFERENCE INTERPOLATION POLYNOMIAL.
C CALLED BY SUBROUTINE FCNSUBZO
C
C EQPAR(K.1) STORE NO. OF DATA POINTS
C EQPAR(K.2) STORE IDEG
C EQPAR(20,3) STORE SLOPE l
c EQPAR(2O.A) STORE SLOPE 2
C

$INSERT COMFILE
$INSERT COMAREA
C ............................................................
lO WRITE(JLUN,900)
CALL CMREAD
IF(NOCMN.EQ.0)THEN
N=5
EQPAR(K,I)IN
GOTO 20
ELSE
DECODE(8O.A.CMN.ERR=899)N
EQPAR(K,I)IN
ENDIF

20 MIN-1
IF(OPT3.EQ.'LINEAR') THEN
IDEGII
ELSE IF(OPT3.EQ.'INTPLO') THEN
IDEGIM
ENDIF
EQPAR(K.2)IIDEG

DO 100 I-I,N
3O WRITE(JLUN,920) 1,1

CALL CMREAD

IF(NOCMN.EQ.0)THEN
XCO(IL.1)-1

ELSE
DECODE(80,A,CMN,ERR=AO)CON
XCO(IL.I)ICON

ENDIF

IF(I.EQ.1)GOTO 50
DUM=XCO(IL.1)-XCO(IL,1-1)
IF(DUM.LE.O.)THEN
WRITE(JLUN.930)
GOTO 10
ELSE
GOTO 50
ENDIF

AO CALL ERRMSG (A)
GOTO 30

50 WRITE(JLUN,9AO) 1,1
CALL CMREAD
IF(NOCMN.EQ.0)THEN
YCO(IL,I)II
ELSE
DECODE(80,A,CMN,ERR=6O)CONI
YCO(IL.I)ICONI

ENDIF
GOTO 100

60 CALL ERRMSG(A)
GOTO 50

100 CONTINUE

C
IF(1L.EQ.2)GDTO 125

C
WRITE(JLUN,950) XCO(IL,1)
CALL GENCON(K,20.3.I)
WRITE(JLUN,960) XCO(IL,N)
CALL GENCON(K,20,A,2)

C

125 CALL DTABLE(1L.N,M)
RETURN

C

899 CALL ERRMSG(A)
GOTO 10

C

C ----- FORMAT -----

C

900 FORMAT(/,‘ENTER NUMBER OF DATA POINTS (5):')

920 FORMAT(/.'ENTER COORDINATE X(',IZ,'), (',12.'):')

930 FORMAT(/,' ERROR] X MUST ARRANGED IN ASCENDING ORDER,DO IT AGAIN')
9A0 FORMAT(/.'ENTER COORDINATE Y(',12,'), (',12,'):')

950 FORMAT(/.'DY/DX FOR X LESS THAN '.E12.A.

+ ' IS Cl')
960 FORMAT(/.'DY/DX FOR x GREATER THAN ',E12.A,
+ ' IS c2')
C
. END
C

C
SUBROUTINE DTABLE(1L.N,M)
C
C ----- CALLED BY SUBROUTINE INTERP
C
COMMON/EQA/ XCO(2.20).YCO(2.20).TABLE(2.20.20)
C ............................................................
IO NMIIN-l
DO 30 I=I.NH1
3O TABLE(IL.I,I)-( YCO(IL,I+l)-YCO(IL,I) )
+ /( XCO(IL,I+I)-XC0(IL.1) )
IF(M.LE.I)GOTO 100
c .
DO 50 J-2,M
DO 50 IIJ,NMI
lSUB-I+l-J
TABLE(IL,I,J)=(TABLE(IL.I.J-1)-TABLE(IL,I-l,J-l))
+ / (XCO(IL,I+l)-XCO(IL,ISUB) )
50 CONTINUE
C
100 CONTINUE
RETURN
END
C
C ............................................................
C
FUNCTION FMIDV(IL,N,IDEG,XARG)
C
C ----- CALLED BY SUBROUTINE NLEVAL, GETVAL
c
COMMON/EQA/ XCO(2,20),YCO(2.20).TABLE(2.2O.20)
C ............................................................
DO 100 1-1,N
IF(I.EQ.N .OR. XARG.LE.XCD(IL,I)) GOTO 200
100 CONTINUE
C
200 MAx-I+IDEG/2
IF(MAX.LE.IDEG) MAXIIDEG+1
1F(MAX.GT.N) MAXIN
C
YEST=TABLE(1L.MAX-I.IDEG)
1F(IDEG.LE.I)GOTO 300
IDEGMI-IDEG-I
DO 250 1=I,IDEGM1
ISUBI-MAX-I
lSUBZ-IDEG-I
250 YEST-YESTA(XARG-XCO(1L,ISUBI)) + TABLE(IL,ISUBl-l,ISUBZ)
C
300 ISUBl-MAX-IDEG
FHIDV=YESTA(XARG—XCO(1L,ISUBI))+YCO(IL,ISUBI)
C

RETURN
END

C 58

C

C---FCNSU822 --------------------------------------------------
SUBROUTINE FCNSUB22(K)

C

C ----- FUNCTION TYPE : RCPSQR

C

C Y . C0 + C1 / ( I+(X/C2)AA2 )

C

$INSERT COMFILE
$INSERT COMAREA

C ............................................................
IF(OPT1.EQ.'LOOK ') GOTO 50
IF(RVS(K).EQ.I)THEN
OPT2I'N'
CALL ERRMSG(6)
RETURN
ENDIF
C

WRITE(JLUN.IOOI)CHIO(K,1),IEQZS(K,1),CHIO(K,2),IEQZS(K,2)
1001 FORMAT(/,5X,AI,12,' I C0 + C1 / ( I + (',AI,12,' / C2)**2 )'./)

CALL GENCON(K.22,1,O)
CALL GENCON(K,22,2.I)
CALL GENCON(K.22.3.2)
CALL ERRMSG(7)
CALL CMREAD
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 50
RETURN
C
50 WRITE(JLUN,IOO3)CHIO(K,1),IEQZS(K,1),EQPAR(K,I),EQPAR(K,2),
+ CHIO(K,2),IEQZS(K,2),EQPAR(K,3)
1003 FORMAT(/,5X,Al,|2,' I',lPE12.A,' + ',1PE12.A,' / ( l + (',

+ Al,12,' /',1PE12.A,')**2 )' )

RETURN

END
C
C---NLENAL ---------------------------------------------------
C

SUBROUTINE NLEVAL
C
C ----- WRITTEN BY : TSUN-YU CHOU, AUGUST 1982
C

$INSERT COMFILE
$INSERT COMAREA
c ............................................................
C
20 WRITE(JLUN,21)
21 FORMAT(/,‘GIVE FUNCTION NUMBER :')
CALL CMREAD
IF(NOCMN.EQ.0)GOTO 20
DECODE(80.A.CMN,ERR-AO)MFCN
IF (MFCN.LT.O) THEN
PRINT A,'ERROR] INCORRECT FUNCTION NUMBER'

GOTO 20

ELSE IF(MFCN.EQ.0)THEN
RETURN

ENDIF

DD 30 JJ-1,NEQS
|F(IEQ2$(JJ,I) .EQ. MFCN) THEN
I=JJ
GOTO 60
ENDIF
3O CONTINUE
AD PRINT A,'ERROR] BAD FUNCTION LABEL.'

GOTO 20
C
60 OPTI='LOOK '
OPT2='S'
C
CALL FCNLOOK(I)
C

100 WRITE(JLUN,101)

lOl FORMAT(/.'GIVE THE VALUE OF INPUT VARIABLE :')
CALL CMREAD
DECODE(80,A.CHN.ERR=120)XXX

CALL GETVAL(XXX,Y,I,IFCNTP(I))
WRITE(JLUN,110)Y
110 FORMAT(/,‘VALUE'OF OUTPUT VARIABLE IS : '.1PE12.A,/)

GO TO 130

C

120 PRINT *,'ERROR, BAD DATA'
GOTO 100

C

130 WRITE(JLUN,1AO)

1A0 FORMAT(/.'EVALUATE ANOTHER FUNCTION 2 (YES):')
CALL CMREAD
IF(NOCMN.EQ.0)GOTO 20
IF(CMN.EQ.'Y' .OR. CMN.EQ.'YES')GOTO 20

C
RETURN
END
C
C---GETVAL ---------------------------------------------------
C
SUBROUTINE GETVAL(XINP,YOUTP,NOEQ,NFCNTP)
C
C ----- WRITTEN BY : TSUN-YU CHOU, AUGUST 1982
C

$INSERT COMFILE
$INSERT COMAREA
C ............................................................
CI
IF(NFCNTP.EQ.I)THEN
YOUTP=EQPAR(NOEQ,1)

C2

C3

CA

C5

C6

C7

+

.1.

+

+

+

+

+

+

RETURN

ELSE IF(NFCNTP.EQ.2)THEN
YOUTPIEQPAR(NOEQ,1) + EQPAR(NOEQ,2) A XINP
RETURN

ELSE IF(NFCNTP.EQ.3)THEN

IF(RVS(NOEQ).EQ.0)THEN

YOUTP-EQPAR(NOEQ.I)ASIN(EQPAR(NOEQ.2)AXINP+EQPAR<NOEO.3))
+EQPAR(NOEQ,A)

ELSE IF(RVS(NOEQ).EQ.I)THEN

YOUTP-EQPAR(NOEQ,I)AASIN(EQPAR(N0EQ,2)AXINP+EQPAR(NOEQ.3))
+EQPAR(NOEQ.A)

ENDIF

RETURN

ELSE IF(NFCNTP.EQ.A)THEN
IF(RVS(NOEQ).EQ.0)THEN
YOUTPIEQPAR(NOEQ,l)ACOS(EQPAR(N0EQ,2)AXINP+EQPAR(NOEQ,3))
+EQPAR(NOEQ,A)
ELSE IF(RVS(NOEQ).EQ.I)THEN
YOUTP-EQPAR(NOEQ,I)AACOS(EQPAR(NOEQ,2)AXINP+EQPAR(NOEQ,3))
+EQPAR(NOEQ,A)
ENDIF
RETURN

ELSE 1F(NFCNTP.EQ.5)THEN

IF(RVS(NOEQ).EQ.0)THEN

YOUTP-EQPAR(NOEQ,I)AASIN(EQPAR(NOEQ,2)AXINP+EQPAR(NOEQ.3))
+EQPAR(NOEQ,A)

ELSE IF(RVS(NOEQ).EQ.I)THEN

YOUTP=EQPAR(N0EQ,I)ASIN(EQPAR(NOEQ,2)AXINP+EQPAR(N0EQ,3))
+EQPAR(NOEQ.A)

ENDIF

RETURN

ELSE 1F(NFCNTP.EQ.6)THEN
IF(RVS(NOEQ).EQ.0)THEN
YOUTP=EQPAR(NOEQ.1)AACOS(EQPAR(NOEQ,2)AXINP+EQPAR(NOEQ,3))
+EQPAR(NOEQ.A)
ELSE IF(RVS(NOEQ).EQ.I)THEN
YOUTPIEQPAR(NOEQ.l)*COS(EQPAR(NOEQ,2)AXINP+EQPAR(NOEQ,3))
+EQPAR(NOEQ.A)
ENDIF
RETURN

ELSE IF(NFCNTP.EQ.7)THEN
IF(RVS(NOEQ).EQ.0)THEN
YOUTP-EQPAR(NOEQ.I) A XINP A ABS(XINP)
ELSE IF(RVS(NOEQ).EQ.I)THEN
YOUTP-EQPAR(NOEQ,I)ASIGN(I,X1NP)ASQRT(ABS(X1NP))
ENDIF
RETURN

C8
ELSE IF(NFCNTP.EQ.8)THEN
IF(RVS(NOEQ).EQ.0)THEN
YOUTP=EQPAR(NOEQ,I)ASIGN(I,X1NP)ASQRT(ABS(XINP))
ELSE IF(RVS(NOEQ).EQ.I)THEN
YOUTP-EQPAR(NOEQ,I) A XINP A ABS(XINP)
ENDIF
RETURN
C9
ELSE IF(NFCNTP.EQ.9)THEN
IF(RVS(NOEQ).EQ.0)THEN
YOUTP-EQPAR(NOEQ.I)AEXP(EQPAR(NOEQ.2)AXINP+EQPAR(NOEQ.3) )
+ +EQPAR(NOEQ,A)
ELSE IF(RVS(NOEQ).EQ.I)THEN
YOUTP=EQPAR(NOEQ,1)AALOG(EQPAR(N0EQ,2)AXINP+EQPAR(NOEQ,3))
+ +£QPAR(NOEQ,A)
ELSE IF(RVS(NOEQ).EQ.I)THEN
ENDIF
RETURN
c10
ELSE IF(NFCNTP.EQ.10)THEN
IF(RVS(NOEQ).EQ.0)THEN
YOUTP=EQPAR(NOEQ.I)AALOG(EQPAR(NOEQ,2)AXINP+EQPAR(NOEQ.3))
+ +EQPAR(NOEQ,A)
ELSE IF(RVS(NOEQ).EQ.I)THEN
YOUTPIEQPAR(NOEQ,1)AEXP(EQPAR(NOEQ,2)AXINP+EQPAR(NOEQ,3) )
+ +EQPAR(NOEQ.A)
ENDIF
RETURN
Cll
ELSE IF(NFCNTP.EQ.II)THEN
IF(RVS(NOEQ).EQ.0)THEN
YOUTP-EQPAR(NOEQ,I)/SIN(EQPAR(NOEQ,2)AXINP+EQPAR(NOEQ.3))
+ +EQPAR(NOEQ,A)
ELSE
YDUTP=EQPAR(NOEQ,I)AASIN(EQPAR(NOEQ.2)/(XINP-EQPAR(NOEQ.3)))
+ +EQPAR(NOEQ,A)
ENDIF
RETURN
C12
ELSE IF(NFCNTP.EQ.12)THEN
IF(RVS(NOEQ).EQ.0)THEN
YOUTP-EQPAR(NOEQ.1)/COS(EQPAR(NOEQ.2)AXINP+EQPAR(NOEQ.3))
+ +EQPAR(NOEQ,A)
ELSE
YOUTPIEQPAR(NOEQ.l)*ACOS(EQPAR(NOEQ.2)/(XINP-EQPAR(NOEQ.3)))
+ +EQPAR(NOEQ.A)
ENDIF
RETURN
C13
ELSE IF(NFCNTP.EQ.13)THEN
IF(RVS(NOEQ).EQ.0)THEN
YOUTP=EQPAR(NOEQ,l)/ASIN(EQPAR(NOEQ,2)AXINP+EQPAR(NOEQ.3))

+ ‘ +EQPAR(NOEQ,A)
ELSE
YOUTP=EQPAR(NOEQ,I)ASIN(EQPAR(NOEQ,2)/(XINP-EQPAR(NOEQ,3)))
+ +EQPAR(NOEQ,A)
ENDIF
RETURN
CIA
ELSE 1F(NFCNTP.EQ.IA)THEN
IF(RVS(NOEQ).EQ.0)THEN
YOUTP=EQPAR(NOEQ,I)/ACOS(EQPAR(NOEQ.2)AXINP+EQPAR(NOEQ,3))
+ +EQPAR(NOEQ,A)
ELSE
YOUTP-EQPAR(NOEQ,I)ACOS(EQPAR(NOEQ,2)/(XINP-EQPAR(NOEQ,3)))
+ +EQPAR(NOEQ,A)
ENDIF
RETURN
C15
ELSE IF(NFCNTP.EQ.I5)THEN
YOUTP=EQPAR(NOEQ,I) +EQPAR(NOEQ.2)AXINP +EQPAR(NOEQ,3)AXINPAA2
+ +EQPAR(NOEQ,A)AXINPAA3 + EQPAR(NOEQ,5)AXINPAAA
RETURN
C16
ELSE IF(NFCNTP.EQ.I6)THEN
1F(X1NP.EQ.O.)THEN
YOUTP-0.
ELSE
YOUTPISIGN(I,XINP) A EQPAR(NOEQ,I)
ENDIF
RETURN
C17
ELSE IF(NFCNTP.EQ.17)THEN
1F(XINP.LE.EQPAR(NOEQ,1))THEN
YOUTP-EQPAR(NOEQ.2)+EQPAR(NOEQ,3)A(XINP-EQPAR(NOEQ,1))
ELSE
YOUTP=EQPAR(NOEQ.2)+EQPAR(NOEQ,A)A(X1NP-EQPAR(NOEQ,1))
ENDIF
RETURN
C18
ELSE IF(NFCNTP.EQ.18)THEN
IF(XINP.LE.EQPAR(NOEQ,I)) THEN
YOUTP-EQPAR(NOEQ,2)
ELSE IF(XINP.GE.EQPAR(NOEQ.3))THEN
YOUTP-EQPAR(NOEQ.A)
ELSE
YOUTP-EQPAR(NOEQ,2) + ( (EQPAR(NOEQ,A)-EQPAR(NOEQ.2))/
+ (EQPAR(NOEQ.3)-EQPAR(NOEQ.I)))A(XINP-EQPAR(NOEQ,1))
ENDIF
RETURN
C19
ELSE IF(NFCNTP.EQ.19)THEN
1F(XINP.LE.EQPAR(N0EQ,1))THEN
YOUTP=EQPAR(NOEQ.2)+EQPAR(NOEQ,5)A(XINP-EQPAR(NOEQ,1))
ELSE IF(XINP.GE.EQPAR(NOEQ.3))THEN

C20

C21

100

110

120

C22

c---

C 63

YOUTP=EQPAR(NOEQ.A)+EQPAR(NOEQ.6)A(XINP-EQPAR(NOEQ.3))
ELSE
YOUTP-EQPAR(NOEQ,2) + ((EQPAR(NOEQ.A)-EQPAR(NOEQ.2)) /
+ (EQPAR(NDEQ.3)-EQPAR(NOEQ.I))) A (XINP-EQPAR(NOEQ,1))
ENDIF
RETURN

ELSE IF(NFCNTP.EQ.20)THEN

IL-I

NN=EQPAR(NOEQ,1)

IDEGIEQPAR(NOEQ,2)

1F((XINP.GE.XCO(IL.I)).AND.(X1NP.LE.XCO(1L,NN)))THEN
YOUTP-FMIDV(IL,NN,1DEG,XINP)

ELSE IF(XINP.LT.XCO(IL,1))THEN
YOUTP=YCO(IL,1)+EQPAR(NOEQ,3)A(XINP-XCO(1L,1))

ELSE IF(XINP.GT.XCO(IL,NN))THEN
Y0UTP=YCO(IL,NN)+EQPAR(NOEQ,A)A(XINP-XCO(IL,NN))

ENDIF

RETURN

ELSE IF(NFCNTP.EQ.21)THEN
ILIZ
NP-EQPAR(NOEQ.I)
IDEGIEQPAR(NOEQ.2)
PRINT A,'ORDER OF POLYNOMIAL CAN BE CHANGED'
PRINT A.'ENTER THE ORDER :'
CALL CMREAD
DECODE(80,A,CMN,ERR-IIO)NI
IDEGINI
IF(IDEG.GE.NP)THEN
WRITE(JLUN,120)NP
GOTO IOO
ENDIF
YOUTP-FMIDV(IL.NP,IDEG,x1NP)
RETURN
PRINT A,'ERROR, BAD DATA.’
GOTO 100
FORMAT(/.'ORDER MUST BE LESS THAN ',13)

ELSE IF(NFCNTP.EQ.22)THEN
YOUTPIEQPAR(NOEQ,I)+ EQPAR(NOEQ,2) / (l+(XINP/EQPAR(NOEQ,3))AAz)
RETURN

RETURN
END

C---COMF|LE

C

INTEGERAA IERRF,MCMAM

CHARACTER*A IALPHA,IELNAM(50),NTYP(9),HEAD(20)
CHARACTER*1 LINE(72),LABEL(72)

CHARACTERAI IEPFQL(20),IXUDL(20),LBL1(1,5)

C
COMMON /BKl/ ILUN,JLUN,IERRF,IPRLST(20).HEAD
COMMON /BK2/ NREAD,LINE,LABEL,INT,REAL,IALPHA,KEY,NCOL,
+ LI,MARGIN,NEXT,IST,ISP
COMMON /BK3/ NEL,NBD,1ELLST(5O),IELNAM,NB1MX(100),
+ IBMX(50,2),NPTR(51),MNEL,MNBD,NTYP
COMMON /BKA/ NMCR,NMCPT,NSTD,MNMCR,MNMCPT,
+ MCNAM(IO),MPTRA(IO),MPTRB(lO),MBIMX(50)
COMMON /BK5/ ICMX(lOO),MELMNT(AO),MBND(AO),M|T(AO),NTEMP(AO),
+ MNSTKM,MNSTKN,MTOP,NTOP,LEVEL,IBLIS(13).ICCNOD(I3),
+ IBPNT(A),NODLIS(A),JPNT,KBOND,INODE,MAXBD.MAXND
COHMDN /BK6/ IBEQ(SD),IBT(5O).NBEX,NFI.NFD,NFL,NFS,NBIN
COMMON /BK7/ PAR(IOO),1PTR(5I),NPAR,MNPAR,IPFLG
COMMON /BK8/ S(150).IPS(150),LENS,HAXS,
+ T(150).IPT(150).LENT.MAXT.
+ WK(lOO),IWK(lOO),LENWK,MAXWK
COMMON /BK9/ FL(IOO).1PFL(IOO).LENFL,MAXFL.
+ FS(100),IPFS(IOO),LENFS,MAXFS,
+ W2(lOO).IW2(100).LENW2,MAXW2
COMMON /BK10/ A(150),1PA(150),LENA,MAXA,
+ B(IOO),IPB(100),LENB.HAXB,
+ E(50),|PE(50).LENE,MAXE.
+ IEPFQL,IBDL(2O),1XUDL.IDXL(20).NREQ,MAXXUD
COMMON /BKll/ RM(225).IPRM(225),LENRM,MAXRM,
+ RN(225).IPRN(225).LENRN,HAXRN,
+ X1N(15).X(15),IPX(15).LENX.MAXX,
+ U(IO),1PU(10),LENU,MAXU,
+ DX(15).IPDX(15),LENDX,
+ lUTYP(lO),UPAR(lO,lO),NARGU(lO),
+ TIN,TF|N,DT,KLIM,DEL,NSAV,MAXRES,RES(101,20)
COMMON /BK12/ T1,T2,DELTA,JX,JDX(5).LBL2(2.5).JDXMAX,
+ YMAX(20),YMIN(20),YSCLO,YSCHI,LBL1
COMMON /BKl3/ WR(10),WI(IO)
COMMON /BKlA/ DATAP(5,IOI)
C
CAAAEND COMFILE
C ......................................................
C
C .......................................................
C-—-COMAREA
C
CHARACTERABO C1NPUT.CMN
CHARACTERA6 CMATCH(3O),CMSETI(3O),1NFSET(20)
CHARACTERA6 FCNSET(3O).OPTI,OPT3
CHARACTERAI CH10(50.A),OPT2.0PTA
INTEGER FPOCMN
C

COMMON/READER/ CINPUT,CMN,FPOCMN,LPOCMN,NOCMN,
CMATCH,IFOUND

C65

COMMON /CMI/ CMSETI,INFSET,FCNSET
COMMON /CTL/ IPNT(IO),OPT1,OPT2,OPT3,OPTA
COMMON /EQI/ NEQS,NFCN(50),NEQTP(50),IEQZS(50,2),
+ CHID.NOID(50)
COMMON /EQ2/ PARDF(30,IO).EQPAR(5O,IO)
COMMON /EQ3/ IFCNTP(5O),RVS(5O)
COMMON /EQA/ XCO(2.2O),YC0(2.2O).TABLE(2.20.20)
C
CAAAEND COMAREA
C ......................................................

 
    

            

IIIIIIIIIIIIIIIIIIIIII

I IIIIIIIIIIIIIIIIIIIIIIIIIIII

3 4

IIIIIIIIII

IIIIIIIIIIILIIII