src/rt/data.c

/* Copyright (c) 1996 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 *  data.c - routines dealing with interpolated data.
 */

#include  "standard.h"

#include  "color.h"

#include  "resolu.h"

#include  "data.h"

                                /* picture memory usage before warning */
#ifndef PSIZWARN
#ifdef BIGMEM
#define PSIZWARN        5000000
#else
#define PSIZWARN        1500000
#endif
#endif

#ifndef TABSIZ
#define TABSIZ          97              /* table size (prime) */
#endif

#define hash(s)         (shash(s)%TABSIZ)


extern char  *getlibpath();             /* library search path */

static DATARRAY  *dtab[TABSIZ];         /* data array list */


DATARRAY *
getdata(dname)                          /* get data array dname */
char  *dname;
{
        char  *dfname;
        FILE  *fp;
        int  asize;
        register int  i, j;
        register DATARRAY  *dp;
                                                /* look for array in list */
        for (dp = dtab[hash(dname)]; dp != NULL; dp = dp->next)
                if (!strcmp(dname, dp->name))
                        return(dp);             /* found! */

        /*
         *      If we haven't loaded the data already, we will look
         *  for it in the directories specified by the library path.
         *
         *      The file has the following format:
         *
         *              N
         *              beg0    end0    n0
         *              beg1    end1    n1
         *              . . .
         *              begN    endN    nN
         *              data, later dimensions changing faster
         *              . . .
         *
         *      For irregularly spaced points, the following can be
         *  substituted for begi endi ni:
         *
         *              0 0 ni p0i p1i .. pni
         */

        if ((dfname = getpath(dname, getlibpath(), R_OK)) == NULL) {
                sprintf(errmsg, "cannot find data file \"%s\"", dname);
                error(USER, errmsg);
        }
        if ((fp = fopen(dfname, "r")) == NULL) {
                sprintf(errmsg, "cannot open data file \"%s\"", dfname);
                error(SYSTEM, errmsg);
        }
                                                        /* get dimensions */
        if (fgetval(fp, 'i', &asize) <= 0)
                goto scanerr;
        if (asize <= 0 | asize > MAXDDIM) {
                sprintf(errmsg, "bad number of dimensions for \"%s\"", dname);
                error(USER, errmsg);
        }
        if ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
                goto memerr;
        dp->name = savestr(dname);
        dp->type = DATATY;
        dp->nd = asize;
        asize = 1;
        for (i = 0; i < dp->nd; i++) {
                if (fgetval(fp, DATATY, &dp->dim[i].org) <= 0)
                        goto scanerr;
                if (fgetval(fp, DATATY, &dp->dim[i].siz) <= 0)
                        goto scanerr;
                if (fgetval(fp, 'i', &dp->dim[i].ne) <= 0)
                        goto scanerr;
                if (dp->dim[i].ne < 2)
                        goto scanerr;
                asize *= dp->dim[i].ne;
                if ((dp->dim[i].siz -= dp->dim[i].org) == 0) {
                        dp->dim[i].p = (DATATYPE *)
                                        malloc(dp->dim[i].ne*sizeof(DATATYPE));
                        if (dp->dim[i].p == NULL)
                                goto memerr;
                        for (j = 0; j < dp->dim[i].ne; j++)
                                if (fgetval(fp, DATATY, &dp->dim[i].p[j]) <= 0)
                                        goto scanerr;
                        for (j = 1; j < dp->dim[i].ne-1; j++)
                                if ((dp->dim[i].p[j-1] < dp->dim[i].p[j]) !=
                                        (dp->dim[i].p[j] < dp->dim[i].p[j+1]))
                                        goto scanerr;
                        dp->dim[i].org = dp->dim[i].p[0];
                        dp->dim[i].siz = dp->dim[i].p[dp->dim[i].ne-1]
                                                - dp->dim[i].p[0];
                } else
                        dp->dim[i].p = NULL;
        }
        if ((dp->arr.d = (DATATYPE *)malloc(asize*sizeof(DATATYPE))) == NULL)
                goto memerr;
        
        for (i = 0; i < asize; i++)
                if (fgetval(fp, DATATY, &dp->arr.d[i]) <= 0)
                        goto scanerr;
        fclose(fp);
        i = hash(dname);
        dp->next = dtab[i];
        return(dtab[i] = dp);

memerr:
        error(SYSTEM, "out of memory in getdata");
scanerr:
        sprintf(errmsg, "%s in data file \"%s\"",
                        feof(fp) ? "unexpected EOF" : "bad format", dfname);
        error(USER, errmsg);
}


static
headaspect(s, iap)                      /* check string for aspect ratio */
char  *s;
double  *iap;
{
        char    fmt[32];

        if (isaspect(s))
                *iap *= aspectval(s);
        else if (formatval(fmt, s) && strcmp(fmt, COLRFMT))
                *iap = 0.0;
        return(0);
}


DATARRAY *
getpict(pname)                          /* get picture pname */
char  *pname;
{
        double  inpaspect;
        char  *pfname;
        FILE  *fp;
        COLR  *scanin;
        int  sl, ns;
        RESOLU  inpres;
        FLOAT  loc[2];
        int  y;
        register int  x, i;
        register DATARRAY  *pp;
                                                /* look for array in list */
        for (pp = dtab[hash(pname)]; pp != NULL; pp = pp->next)
                if (!strcmp(pname, pp->name))
                        return(pp);             /* found! */

        if ((pfname = getpath(pname, getlibpath(), R_OK)) == NULL) {
                sprintf(errmsg, "cannot find picture file \"%s\"", pname);
                error(USER, errmsg);
        }
        if ((pp = (DATARRAY *)malloc(3*sizeof(DATARRAY))) == NULL)
                goto memerr;

        pp[0].name = savestr(pname);

        if ((fp = fopen(pfname, "r")) == NULL) {
                sprintf(errmsg, "cannot open picture file \"%s\"", pfname);
                error(SYSTEM, errmsg);
        }
#ifdef MSDOS
        setmode(fileno(fp), O_BINARY);
#endif
                                                /* get dimensions */
        inpaspect = 1.0;
        getheader(fp, headaspect, &inpaspect);
        if (inpaspect <= FTINY || !fgetsresolu(&inpres, fp))
                goto readerr;
        pp[0].nd = 2;
        pp[0].dim[0].ne = inpres.yr;
        pp[0].dim[1].ne = inpres.xr;
        pp[0].dim[0].org =
        pp[0].dim[1].org = 0.0;
        if (inpres.xr <= inpres.yr*inpaspect) {
                pp[0].dim[0].siz = inpaspect *
                                        (double)inpres.yr/inpres.xr;
                pp[0].dim[1].siz = 1.0;
        } else {
                pp[0].dim[0].siz = 1.0;
                pp[0].dim[1].siz = (double)inpres.xr/inpres.yr /
                                        inpaspect;
        }
        pp[0].dim[0].p = pp[0].dim[1].p = NULL;
        sl = scanlen(&inpres);                          /* allocate array */
        ns = numscans(&inpres);
        i = ns*sl*sizeof(COLR);
#if PSIZWARN
        if (i > PSIZWARN) {                             /* memory warning */
                sprintf(errmsg, "picture file \"%s\" using %d bytes of memory",
                                pname, i);
                error(WARNING, errmsg);
        }
#endif
        if ((pp[0].arr.c = (COLR *)malloc(i)) == NULL)
                goto memerr;
                                                        /* load picture */
        if ((scanin = (COLR *)malloc(sl*sizeof(COLR))) == NULL)
                goto memerr;
        for (y = 0; y < ns; y++) {
                if (freadcolrs(scanin, sl, fp) < 0)
                        goto readerr;
                for (x = 0; x < sl; x++) {
                        pix2loc(loc, &inpres, x, y);
                        i = (int)(loc[1]*inpres.yr)*inpres.xr +
                                        (int)(loc[0]*inpres.xr);
                        copycolr(pp[0].arr.c[i], scanin[x]);
                }
        }
        free((char *)scanin);
        fclose(fp);
        i = hash(pname);
        pp[0].next = dtab[i];           /* link into picture list */
        copystruct(&pp[1], &pp[0]);
        copystruct(&pp[2], &pp[0]);
        pp[0].type = RED;               /* differentiate RGB records */
        pp[1].type = GRN;
        pp[2].type = BLU;
        return(dtab[i] = pp);

memerr:
        error(SYSTEM, "out of memory in getpict");
readerr:
        sprintf(errmsg, "bad picture file \"%s\"", pfname);
        error(USER, errmsg);
}


freedata(dname)                 /* free memory associated with dname */
char  *dname;
{
        DATARRAY  head;
        int  hval, nents;
        register DATARRAY  *dp, *dpl;
        register int  i;

        if (dname == NULL) {                    /* free all if NULL */
                hval = 0; nents = TABSIZ;
        } else {
                hval = hash(dname); nents = 1;
        }
        while (nents--) {
                head.next = dtab[hval];
                dpl = &head;
                while ((dp = dpl->next) != NULL)
                        if (dname == NULL || !strcmp(dname, dp->name)) {
                                dpl->next = dp->next;
                                if (dp->type == DATATY)
                                        free((char *)dp->arr.d);
                                else
                                        free((char *)dp->arr.c);
                                for (i = 0; i < dp->nd; i++)
                                        if (dp->dim[i].p != NULL)
                                                free((char *)dp->dim[i].p);
                                freestr(dp->name);
                                free((char *)dp);
                        } else
                                dpl = dp;
                dtab[hval++] = head.next;
        }
}


double
datavalue(dp, pt)               /* interpolate data value at a point */
register DATARRAY  *dp;
double  *pt;
{
        DATARRAY  sd;
        int  asize;
        int  lower, upper;
        register int  i;
        double  x, y0, y1;
                                        /* set up dimensions for recursion */
        if (dp->nd > 1) {
                sd.name = dp->name;
                sd.type = dp->type;
                sd.nd = dp->nd - 1;
                asize = 1;
                for (i = 0; i < sd.nd; i++) {
                        sd.dim[i].org = dp->dim[i+1].org;
                        sd.dim[i].siz = dp->dim[i+1].siz;
                        sd.dim[i].p = dp->dim[i+1].p;
                        asize *= sd.dim[i].ne = dp->dim[i+1].ne;
                }
        }
                                        /* get independent variable */
        if (dp->dim[0].p == NULL) {             /* evenly spaced points */
                x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
                x *= (double)(dp->dim[0].ne - 1);
                i = x;
                if (i < 0)
                        i = 0;
                else if (i > dp->dim[0].ne - 2)
                        i = dp->dim[0].ne - 2;
        } else {                                /* unevenly spaced points */
                if (dp->dim[0].siz > 0.0) {
                        lower = 0;
                        upper = dp->dim[0].ne;
                } else {
                        lower = dp->dim[0].ne;
                        upper = 0;
                }
                do {
                        i = (lower + upper) >> 1;
                        if (pt[0] >= dp->dim[0].p[i])
                                lower = i;
                        else
                                upper = i;
                } while (i != (lower + upper) >> 1);
                if (i > dp->dim[0].ne - 2)
                        i = dp->dim[0].ne - 2;
                x = i + (pt[0] - dp->dim[0].p[i]) /
                                (dp->dim[0].p[i+1] - dp->dim[0].p[i]);
        }
                                        /* get dependent variable */
        if (dp->nd > 1) {
                if (dp->type == DATATY) {
                        sd.arr.d = dp->arr.d + i*asize;
                        y0 = datavalue(&sd, pt+1);
                        sd.arr.d = dp->arr.d + (i+1)*asize;
                        y1 = datavalue(&sd, pt+1);
                } else {
                        sd.arr.c = dp->arr.c + i*asize;
                        y0 = datavalue(&sd, pt+1);
                        sd.arr.c = dp->arr.c + (i+1)*asize;
                        y1 = datavalue(&sd, pt+1);
                }
        } else {
                if (dp->type == DATATY) {
                        y0 = dp->arr.d[i];
                        y1 = dp->arr.d[i+1];
                } else {
                        y0 = colrval(dp->arr.c[i],dp->type);
                        y1 = colrval(dp->arr.c[i+1],dp->type);
                }
        }
        /*
         * Extrapolate as far as one division, then
         * taper off harmonically to zero.
         */
        if (x > i+2)
                return( (2*y1-y0)/(x-(i-1)) );

        if (x < i-1)
                return( (2*y0-y1)/(i-x) );

        return( y0*((i+1)-x) + y1*(x-i) );
}
Revision:	2.14
Committed:	Thu Nov 5 17:11:56 1998 UTC (25 years, 5 months ago) by gwlarson
Content type:	text/plain
Branch:	MAIN
Changes since 2.13:	+6 -1 lines
Log Message:	changed getheader() to pay attention to function return value
#	User	Rev	Content
1	greg	2.12	/* Copyright (c) 1996 Regents of the University of California */
2	greg	1.1
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6
7			/*
8			* data.c - routines dealing with interpolated data.
9			*/
10
11			#include "standard.h"
12
13			#include "color.h"
14
15	greg	1.11	#include "resolu.h"
16
17	greg	1.1	#include "data.h"
18
19	greg	2.7	/* picture memory usage before warning */
20			#ifndef PSIZWARN
21			#ifdef BIGMEM
22	gregl	2.13	#define PSIZWARN 5000000
23	greg	2.7	#else
24	gregl	2.13	#define PSIZWARN 1500000
25	greg	2.7	#endif
26			#endif
27	greg	1.1
28	greg	2.7	#ifndef TABSIZ
29	greg	2.6	#define TABSIZ 97 /* table size (prime) */
30	greg	2.7	#endif
31	greg	2.6
32			#define hash(s) (shash(s)%TABSIZ)
33
34
35	greg	2.8	extern char getlibpath(); / library search path */
36	greg	1.1
37	greg	2.6	static DATARRAY dtab[TABSIZ]; / data array list */
38	greg	1.1
39
40			DATARRAY *
41			getdata(dname) /* get data array dname */
42			char *dname;
43			{
44			char *dfname;
45			FILE *fp;
46			int asize;
47	greg	1.4	register int i, j;
48	greg	1.1	register DATARRAY *dp;
49			/* look for array in list */
50	greg	2.6	for (dp = dtab[hash(dname)]; dp != NULL; dp = dp->next)
51	greg	1.1	if (!strcmp(dname, dp->name))
52			return(dp); /* found! */
53
54			/*
55			* If we haven't loaded the data already, we will look
56	greg	2.5	* for it in the directories specified by the library path.
57	greg	1.1	*
58			* The file has the following format:
59			*
60	greg	1.4	* N
61	greg	1.1	* beg0 end0 n0
62			* beg1 end1 n1
63			* . . .
64	greg	1.4	* begN endN nN
65	greg	1.1	* data, later dimensions changing faster
66			* . . .
67			*
68	greg	1.4	* For irregularly spaced points, the following can be
69			* substituted for begi endi ni:
70			*
71	greg	1.6	* 0 0 ni p0i p1i .. pni
72	greg	1.1	*/
73
74	greg	2.8	if ((dfname = getpath(dname, getlibpath(), R_OK)) == NULL) {
75	greg	1.1	sprintf(errmsg, "cannot find data file \"%s\"", dname);
76			error(USER, errmsg);
77			}
78			if ((fp = fopen(dfname, "r")) == NULL) {
79			sprintf(errmsg, "cannot open data file \"%s\"", dfname);
80			error(SYSTEM, errmsg);
81			}
82			/* get dimensions */
83	greg	2.12	if (fgetval(fp, 'i', &asize) <= 0)
84	greg	1.1	goto scanerr;
85	greg	2.12	if (asize <= 0 \| asize > MAXDDIM) {
86	greg	1.1	sprintf(errmsg, "bad number of dimensions for \"%s\"", dname);
87			error(USER, errmsg);
88			}
89	greg	2.12	if ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
90			goto memerr;
91			dp->name = savestr(dname);
92			dp->type = DATATY;
93			dp->nd = asize;
94	greg	1.1	asize = 1;
95			for (i = 0; i < dp->nd; i++) {
96	greg	2.12	if (fgetval(fp, DATATY, &dp->dim[i].org) <= 0)
97	greg	1.6	goto scanerr;
98	greg	2.12	if (fgetval(fp, DATATY, &dp->dim[i].siz) <= 0)
99	greg	1.10	goto scanerr;
100	greg	2.10	if (fgetval(fp, 'i', &dp->dim[i].ne) <= 0)
101	greg	1.10	goto scanerr;
102	greg	1.6	if (dp->dim[i].ne < 2)
103			goto scanerr;
104			asize *= dp->dim[i].ne;
105			if ((dp->dim[i].siz -= dp->dim[i].org) == 0) {
106	greg	2.12	dp->dim[i].p = (DATATYPE *)
107			malloc(dp->dim[i].ne*sizeof(DATATYPE));
108	greg	1.4	if (dp->dim[i].p == NULL)
109			goto memerr;
110	greg	2.10	for (j = 0; j < dp->dim[i].ne; j++)
111	greg	2.12	if (fgetval(fp, DATATY, &dp->dim[i].p[j]) <= 0)
112	greg	1.4	goto scanerr;
113			for (j = 1; j < dp->dim[i].ne-1; j++)
114			if ((dp->dim[i].p[j-1] < dp->dim[i].p[j]) !=
115			(dp->dim[i].p[j] < dp->dim[i].p[j+1]))
116			goto scanerr;
117			dp->dim[i].org = dp->dim[i].p[0];
118			dp->dim[i].siz = dp->dim[i].p[dp->dim[i].ne-1]
119			- dp->dim[i].p[0];
120			} else
121	greg	1.6	dp->dim[i].p = NULL;
122	greg	1.1	}
123	greg	2.12	if ((dp->arr.d = (DATATYPE )malloc(asizesizeof(DATATYPE))) == NULL)
124	greg	1.1	goto memerr;
125
126	greg	2.10	for (i = 0; i < asize; i++)
127	greg	2.12	if (fgetval(fp, DATATY, &dp->arr.d[i]) <= 0)
128	greg	1.1	goto scanerr;
129			fclose(fp);
130	greg	2.6	i = hash(dname);
131			dp->next = dtab[i];
132			return(dtab[i] = dp);
133	greg	1.1
134			memerr:
135			error(SYSTEM, "out of memory in getdata");
136			scanerr:
137			sprintf(errmsg, "%s in data file \"%s\"",
138			feof(fp) ? "unexpected EOF" : "bad format", dfname);
139			error(USER, errmsg);
140			}
141
142
143	greg	1.5	static
144	greg	2.4	headaspect(s, iap) /* check string for aspect ratio */
145	greg	1.5	char *s;
146	greg	2.4	double *iap;
147	greg	1.5	{
148	gwlarson	2.14	char fmt[32];
149
150	greg	1.5	if (isaspect(s))
151	greg	2.4	iap = aspectval(s);
152	gwlarson	2.14	else if (formatval(fmt, s) && strcmp(fmt, COLRFMT))
153			*iap = 0.0;
154			return(0);
155	greg	1.5	}
156
157
158	greg	1.1	DATARRAY *
159			getpict(pname) /* get picture pname */
160			char *pname;
161			{
162	greg	2.4	double inpaspect;
163	greg	1.1	char *pfname;
164			FILE *fp;
165	greg	2.12	COLR *scanin;
166	greg	1.11	int sl, ns;
167	greg	2.3	RESOLU inpres;
168	greg	1.11	FLOAT loc[2];
169			int y;
170			register int x, i;
171	greg	1.1	register DATARRAY *pp;
172			/* look for array in list */
173	greg	2.12	for (pp = dtab[hash(pname)]; pp != NULL; pp = pp->next)
174	greg	1.1	if (!strcmp(pname, pp->name))
175			return(pp); /* found! */
176
177	greg	2.8	if ((pfname = getpath(pname, getlibpath(), R_OK)) == NULL) {
178	greg	1.1	sprintf(errmsg, "cannot find picture file \"%s\"", pname);
179			error(USER, errmsg);
180			}
181	greg	2.12	if ((pp = (DATARRAY )malloc(3sizeof(DATARRAY))) == NULL)
182	greg	1.1	goto memerr;
183
184	greg	2.12	pp[0].name = savestr(pname);
185	greg	1.1
186			if ((fp = fopen(pfname, "r")) == NULL) {
187			sprintf(errmsg, "cannot open picture file \"%s\"", pfname);
188			error(SYSTEM, errmsg);
189			}
190	greg	2.3	#ifdef MSDOS
191			setmode(fileno(fp), O_BINARY);
192			#endif
193	greg	1.1	/* get dimensions */
194	greg	1.5	inpaspect = 1.0;
195	greg	2.4	getheader(fp, headaspect, &inpaspect);
196	gwlarson	2.14	if (inpaspect <= FTINY \|\| !fgetsresolu(&inpres, fp))
197	greg	1.1	goto readerr;
198	greg	2.12	pp[0].nd = 2;
199			pp[0].dim[0].ne = inpres.yr;
200			pp[0].dim[1].ne = inpres.xr;
201			pp[0].dim[0].org =
202			pp[0].dim[1].org = 0.0;
203			if (inpres.xr <= inpres.yr*inpaspect) {
204			pp[0].dim[0].siz = inpaspect *
205			(double)inpres.yr/inpres.xr;
206			pp[0].dim[1].siz = 1.0;
207			} else {
208			pp[0].dim[0].siz = 1.0;
209			pp[0].dim[1].siz = (double)inpres.xr/inpres.yr /
210			inpaspect;
211			}
212			pp[0].dim[0].p = pp[0].dim[1].p = NULL;
213			sl = scanlen(&inpres); /* allocate array */
214			ns = numscans(&inpres);
215			i = nsslsizeof(COLR);
216	greg	2.7	#if PSIZWARN
217	greg	2.12	if (i > PSIZWARN) { /* memory warning */
218	greg	2.7	sprintf(errmsg, "picture file \"%s\" using %d bytes of memory",
219			pname, i);
220			error(WARNING, errmsg);
221			}
222			#endif
223	greg	2.12	if ((pp[0].arr.c = (COLR *)malloc(i)) == NULL)
224			goto memerr;
225	greg	1.1	/* load picture */
226	greg	2.12	if ((scanin = (COLR )malloc(slsizeof(COLR))) == NULL)
227	greg	1.1	goto memerr;
228	greg	1.11	for (y = 0; y < ns; y++) {
229	greg	2.12	if (freadcolrs(scanin, sl, fp) < 0)
230	greg	1.1	goto readerr;
231	greg	1.11	for (x = 0; x < sl; x++) {
232			pix2loc(loc, &inpres, x, y);
233			i = (int)(loc[1]inpres.yr)inpres.xr +
234			(int)(loc[0]*inpres.xr);
235	greg	2.12	copycolr(pp[0].arr.c[i], scanin[x]);
236	greg	1.11	}
237	greg	1.1	}
238			free((char *)scanin);
239			fclose(fp);
240	greg	2.6	i = hash(pname);
241	greg	2.12	pp[0].next = dtab[i]; /* link into picture list */
242			copystruct(&pp[1], &pp[0]);
243			copystruct(&pp[2], &pp[0]);
244			pp[0].type = RED; /* differentiate RGB records */
245			pp[1].type = GRN;
246			pp[2].type = BLU;
247			return(dtab[i] = pp);
248	greg	1.1
249			memerr:
250			error(SYSTEM, "out of memory in getpict");
251			readerr:
252			sprintf(errmsg, "bad picture file \"%s\"", pfname);
253			error(USER, errmsg);
254			}
255
256
257			freedata(dname) /* free memory associated with dname */
258			char *dname;
259			{
260	greg	2.5	DATARRAY head;
261	greg	2.6	int hval, nents;
262	greg	1.1	register DATARRAY dp, dpl;
263	greg	1.4	register int i;
264	greg	1.1
265	greg	2.6	if (dname == NULL) { /* free all if NULL */
266			hval = 0; nents = TABSIZ;
267			} else {
268			hval = hash(dname); nents = 1;
269			}
270			while (nents--) {
271			head.next = dtab[hval];
272			dpl = &head;
273			while ((dp = dpl->next) != NULL)
274			if (dname == NULL \|\| !strcmp(dname, dp->name)) {
275			dpl->next = dp->next;
276	greg	2.12	if (dp->type == DATATY)
277			free((char *)dp->arr.d);
278			else
279			free((char *)dp->arr.c);
280	greg	2.6	for (i = 0; i < dp->nd; i++)
281			if (dp->dim[i].p != NULL)
282			free((char *)dp->dim[i].p);
283			freestr(dp->name);
284			free((char *)dp);
285			} else
286			dpl = dp;
287			dtab[hval++] = head.next;
288			}
289	greg	1.1	}
290
291
292			double
293			datavalue(dp, pt) /* interpolate data value at a point */
294			register DATARRAY *dp;
295	greg	2.3	double *pt;
296	greg	1.1	{
297			DATARRAY sd;
298			int asize;
299	greg	1.6	int lower, upper;
300	greg	1.1	register int i;
301	greg	2.9	double x, y0, y1;
302	greg	1.4	/* set up dimensions for recursion */
303	greg	2.12	if (dp->nd > 1) {
304			sd.name = dp->name;
305			sd.type = dp->type;
306			sd.nd = dp->nd - 1;
307			asize = 1;
308			for (i = 0; i < sd.nd; i++) {
309			sd.dim[i].org = dp->dim[i+1].org;
310			sd.dim[i].siz = dp->dim[i+1].siz;
311			sd.dim[i].p = dp->dim[i+1].p;
312			asize *= sd.dim[i].ne = dp->dim[i+1].ne;
313			}
314	greg	1.1	}
315	greg	1.4	/* get independent variable */
316			if (dp->dim[0].p == NULL) { /* evenly spaced points */
317			x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
318	greg	2.9	x *= (double)(dp->dim[0].ne - 1);
319	greg	1.4	i = x;
320			if (i < 0)
321			i = 0;
322			else if (i > dp->dim[0].ne - 2)
323			i = dp->dim[0].ne - 2;
324			} else { /* unevenly spaced points */
325	greg	1.6	if (dp->dim[0].siz > 0.0) {
326			lower = 0;
327			upper = dp->dim[0].ne;
328			} else {
329			lower = dp->dim[0].ne;
330			upper = 0;
331			}
332			do {
333			i = (lower + upper) >> 1;
334			if (pt[0] >= dp->dim[0].p[i])
335			lower = i;
336	greg	1.8	else
337	greg	1.6	upper = i;
338			} while (i != (lower + upper) >> 1);
339	greg	1.8	if (i > dp->dim[0].ne - 2)
340	greg	1.4	i = dp->dim[0].ne - 2;
341			x = i + (pt[0] - dp->dim[0].p[i]) /
342			(dp->dim[0].p[i+1] - dp->dim[0].p[i]);
343			}
344			/* get dependent variable */
345	greg	2.12	if (dp->nd > 1) {
346			if (dp->type == DATATY) {
347			sd.arr.d = dp->arr.d + i*asize;
348			y0 = datavalue(&sd, pt+1);
349			sd.arr.d = dp->arr.d + (i+1)*asize;
350			y1 = datavalue(&sd, pt+1);
351			} else {
352			sd.arr.c = dp->arr.c + i*asize;
353			y0 = datavalue(&sd, pt+1);
354			sd.arr.c = dp->arr.c + (i+1)*asize;
355			y1 = datavalue(&sd, pt+1);
356			}
357	greg	1.1	} else {
358	greg	2.12	if (dp->type == DATATY) {
359			y0 = dp->arr.d[i];
360			y1 = dp->arr.d[i+1];
361			} else {
362			y0 = colrval(dp->arr.c[i],dp->type);
363			y1 = colrval(dp->arr.c[i+1],dp->type);
364			}
365	greg	1.1	}
366			/*
367			* Extrapolate as far as one division, then
368			* taper off harmonically to zero.
369			*/
370			if (x > i+2)
371	greg	2.9	return( (2*y1-y0)/(x-(i-1)) );
372	greg	1.1
373	greg	2.9	if (x < i-1)
374			return( (2*y0-y1)/(i-x) );
375
376			return( y0((i+1)-x) + y1(x-i) );
377	greg	1.1	}