Initial commit

1 files changed, 132 insertions, 0 deletions

diff --git a/pkg/proto/interp.x b/pkg/proto/interp.x
new file mode 100644
index 00000000..026c1c66
--- /dev/null
+++ b/pkg/proto/interp.x
@@ -0,0 +1,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