diff options
Diffstat (limited to 'pkg/utilities/nttools/tintegrate')
| -rw-r--r-- | pkg/utilities/nttools/tintegrate/mkpkg | 11 | ||||
| -rw-r--r-- | pkg/utilities/nttools/tintegrate/tintegrate.x | 155 |
2 files changed, 166 insertions, 0 deletions
diff --git a/pkg/utilities/nttools/tintegrate/mkpkg b/pkg/utilities/nttools/tintegrate/mkpkg new file mode 100644 index 00000000..cfe87a34 --- /dev/null +++ b/pkg/utilities/nttools/tintegrate/mkpkg @@ -0,0 +1,11 @@ +# Update the tintegrate application code in the ttools package library +# Author: GIARETTA, 7-DEC-1987 + +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +libpkg.a: + tintegrate.x <error.h> <ctype.h> <tbset.h> + ; diff --git a/pkg/utilities/nttools/tintegrate/tintegrate.x b/pkg/utilities/nttools/tintegrate/tintegrate.x new file mode 100644 index 00000000..3f01fd67 --- /dev/null +++ b/pkg/utilities/nttools/tintegrate/tintegrate.x @@ -0,0 +1,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 |