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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
|
include <error.h>
include <fset.h> # to check whether input is redirected
include <ctype.h>
include <tbset.h>
# T_INTEGRATE -- integrate one column of a table wrt another
# using simple trapezoid rule, ignoring INDEF values
#
# D. Giaretta, 01-Aug-1987 Original SPP version
# Phil Hodge, 7-Sep-1988 Change parameter name for table.
# Phil Hodge, 26-Jan-1996 Add option to just sum the values.
# Phil hodge, 8-Jun-1999 Set input to STDIN if redirected;
# open the table READ_ONLY, not READ_WRITE.
procedure t_tintegrate()
char inname[SZ_FNAME] # input table name
char col1[SZ_COLNAME] # integrand
char col2[SZ_COLNAME] # independent var.
#--
pointer sp
pointer pt1, pt2, nullpt1, nullpt2
pointer tp
pointer colptr1, colptr2
int i
long numrows
long numused
long lastgood, firstgood # zero indexed
double integ
int fstati()
int tbpsta()
pointer tbtopn()
bool isblank()
errchk tbtopn, clgstr, tbpsta, clputr, clputl
begin
call smark(sp)
integ = 0.0d0
numused = 0
lastgood = -1
firstgood = -1
if (fstati (STDIN, F_REDIR) == YES)
call strcpy ("STDIN", inname, SZ_FNAME)
else
call clgstr ("table", inname, SZ_FNAME)
tp = tbtopn( inname, READ_ONLY, 0)
numrows = tbpsta(tp, TBL_NROWS)
call clgstr( "integrand", col1, SZ_COLNAME)
call clgstr( "independent", col2, SZ_COLNAME)
call tbcfnd( tp, col1, colptr1, 1)
if (colptr1 == NULL)
call error (0, "integrand not found in table")
if (isblank (col2)) {
colptr2 = NULL
} else {
call tbcfnd( tp, col2, colptr2, 1)
if (colptr2 == NULL)
call error (0, "independent variable not found in table")
}
# Get dependent variable values.
call salloc( pt1, numrows, TY_DOUBLE)
call salloc( nullpt1, numrows, TY_BOOL)
call tbcgtd(tp, colptr1, Memd[pt1], Memb[nullpt1], 1, numrows)
# Get independent variable values.
if (colptr2 != NULL) {
call salloc( pt2, numrows, TY_DOUBLE)
call salloc( nullpt2, numrows, TY_BOOL)
call tbcgtd(tp, colptr2, Memd[pt2], Memb[nullpt2], 1, numrows)
}
# Find first non-INDEF row.
if (colptr2 == NULL) {
do i = 0, numrows-1 {
if (!Memb[nullpt1+i]) {
firstgood = i
break
}
}
} else {
do i = 0, numrows-1 {
if ( !Memb[nullpt1+i] && !Memb[nullpt2+i] ) {
firstgood = i
break
}
}
}
if (firstgood == -1) {
call tbtclo (tp)
call error (1, "no data in table")
}
lastgood = firstgood
# apply simple trapezoid rule -
# ignore INDEF values, in other words linearly interpolate
# between adjacent good values
# note also that the independent must be non-decreasing
if (colptr2 == NULL) {
# No independent variable; just sum the values.
numused = 0
do i= firstgood, numrows-1 {
if ( !Memb[nullpt1+i] ) {
integ = integ + Memd[pt1+i]
lastgood = i
numused = numused + 1
}
}
} else {
# Integrate with respect to an independent variable.
numused = 1
do i= firstgood+1, numrows-1 {
if ( !Memb[nullpt1+i] && !Memb[nullpt2+i] ) {
if ( Memd[pt2+i] < Memd[pt2+lastgood] )
call error(0, "independent variable not increasing")
integ = integ +
0.5*(Memd[pt1+i] + Memd[pt1+lastgood]) *
(Memd[pt2+i] - Memd[pt2+lastgood])
lastgood = i
numused = numused + 1
}
}
}
# output integral both as parameter and to STDOUT
# also record the number of good points used
if (numused < 2) {
call printf (" integral = INDEF, at least 2 good points required")
call clputd ("integral", INDEFD)
} else {
call printf (" integral= %g using %d points\n")
call pargd (integ)
call pargi (numused)
call clputd ("integral", integ)
}
call clputl ("ptsused", numused)
call tbtclo (tp)
call sfree (sp)
end
|
