1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
include <fset.h>
define SZ_TABLE 4096
define LINEAR 1
define SPLINE 2
# T_INTERP -- Interpolate for values in a table
#
# A table of x,y pairs contained in a file is used to
# find interpolated values, y, for any other given independent
# variable, x. Extrapolation is performed if necessary.
#
# A series of values may be generated to generate a fine grid
# through a coarse sampling for purposes of plotting. This is
# done by setting the hidden parameter curve_gen to yes.
# The starting point, ending point, and sampling interval
# are also needed in this case (x1, x2, dx).
#
# If only a small number of values are needed to be interpolated
# from the table, the user may enter a number of x's from either
# a file or STDIN.
procedure t_interp()
double x, y, x1, x2, dx
pointer xtab, ytab
int npts, ierr, tbsize
int filelist, tbl, in, imode
char fname[SZ_FNAME], tbl_file[SZ_FNAME], mode[SZ_FNAME]
bool gen
int clpopni(), clgfil(), open(), fscan(), strncmp(), nscan()
real clgetr()
bool clgetb()
begin
# Initialize interpolator
call intrp0 (1)
# File containing x,y pairs in a table
call clgstr ("tbl_file", tbl_file, SZ_FNAME)
# Open table file and read as many points as possible
tbl = open (tbl_file, READ_ONLY, TEXT_FILE)
npts = 0
# Allocate the initial arrays.
call calloc (xtab, SZ_TABLE, TY_DOUBLE)
call calloc (ytab, SZ_TABLE, TY_DOUBLE)
tbsize = SZ_TABLE
while (fscan(tbl) != EOF) {
npts = npts + 1
if (npts > tbsize) {
call realloc (xtab, (tbsize+SZ_TABLE), TY_DOUBLE)
call realloc (ytab, (tbsize+SZ_TABLE), TY_DOUBLE)
tbsize = tbsize + SZ_TABLE
}
call gargd (Memd[xtab+npts-1])
call gargd (Memd[ytab+npts-1])
if (nscan() < 2) {
call eprintf ("Error reading x,y pairs\n")
npts = npts - 1
}
}
call close (tbl)
if (npts < 1)
call error (1, "Table has no entries.")
# Linear or spline interpolator?
call clgstr ("int_mode", mode, SZ_FNAME)
if (strncmp (mode, "linear", 6) == 0)
imode = LINEAR
else
imode = SPLINE
# Generate a curve?
gen = clgetb ("curve_gen")
if (gen) {
x1 = double(clgetr ("x1"))
x2 = double(clgetr ("x2"))
dx = double(clgetr ("dx"))
# Verify that dx will not cause an infinite loop
if (dx == 0.0 || dx * (x2-x1) < 0.0)
call error (1, "Interval paramater dx implies infinite loop.")
for (x=x1; x <= x2; x = x+dx) {
call intrp (1, Memd[xtab], Memd[ytab], npts, x, y, ierr)
call printf ("%12.5g %12.5g\n")
call pargd (x)
call pargd (y)
}
# No, just one point at a time
} else {
# Open input list
filelist = clpopni ("input")
while (clgfil (filelist, fname, SZ_FNAME) != EOF) {
call fseti (STDOUT, F_FLUSHNL, YES)
in = open (fname, READ_ONLY, TEXT_FILE)
# Process input requests
while (fscan(in) != EOF) {
call gargd (x)
if (imode == LINEAR)
call lintrp (1, Memd[xtab], Memd[ytab], npts, x,y, ierr)
else
call intrp (1, Memd[xtab], Memd[ytab], npts, x,y, ierr)
call printf ("%12.5g %12.5g\n")
call pargd (x)
call pargd (y)
}
call close (in)
}
call clpcls (filelist)
}
# Free the pointers.
call mfree (xtab, TY_DOUBLE)
call mfree (ytab, TY_DOUBLE)
end
|
