src/rt/data.c

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

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

/*
 *  data.c - routines dealing with interpolated data.
 *
 *     6/4/86
 */

#include  "standard.h"

#include  "color.h"

#include  "resolu.h"

#include  "data.h"

                                /* picture memory usage before warning */
#ifndef PSIZWARN
#ifdef BIGMEM
#define PSIZWARN        3000000
#else
#define PSIZWARN        1000000
#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 */

static DATARRAY  *ptab[TABSIZ];         /* picture 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 ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
                goto memerr;

        dp->name = savestr(dname);

        if ((fp = fopen(dfname, "r")) == NULL) {
                sprintf(errmsg, "cannot open data file \"%s\"", dfname);
                error(SYSTEM, errmsg);
        }
                                                        /* get dimensions */
        if (fgetval(fp, 'i', &dp->nd) <= 0)
                goto scanerr;
        if (dp->nd <= 0 || dp->nd > MAXDDIM) {
                sprintf(errmsg, "bad number of dimensions for \"%s\"", dname);
                error(USER, errmsg);
        }
        asize = 1;
        for (i = 0; i < dp->nd; i++) {
                if (fgetval(fp, 'd', &dp->dim[i].org) <= 0)
                        goto scanerr;
                if (fgetval(fp, 'd', &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 = (double *)malloc(dp->dim[i].ne*sizeof(double));
                        if (dp->dim[i].p == NULL)
                                goto memerr;
                        for (j = 0; j < dp->dim[i].ne; j++)
                                if (fgetval(fp, 'd', &dp->dim[i].p[i]) <= 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 = (DATATYPE *)malloc(asize*sizeof(DATATYPE))) == NULL)
                goto memerr;
        
        for (i = 0; i < asize; i++)
                if (fgetval(fp, DATATY, &dp->arr[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;
{
        if (isaspect(s))
                *iap *= aspectval(s);
}


DATARRAY *
getpict(pname)                          /* get picture pname */
char  *pname;
{
        double  inpaspect;
        char  *pfname;
        FILE  *fp;
        COLOR  *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 = ptab[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 *)calloc(3, sizeof(DATARRAY))) == NULL)
                goto memerr;

        pp[0].name = 
        pp[1].name =
        pp[2].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 (!fgetsresolu(&inpres, fp))
                goto readerr;
#if PSIZWARN
                                                /* check memory usage */
        i = 3*sizeof(DATATYPE)*inpres.xr*inpres.yr;
        if (i > PSIZWARN) {
                sprintf(errmsg, "picture file \"%s\" using %d bytes of memory",
                                pname, i);
                error(WARNING, errmsg);
        }
#endif
        for (i = 0; i < 3; i++) {
                pp[i].nd = 2;
                pp[i].dim[0].ne = inpres.yr;
                pp[i].dim[1].ne = inpres.xr;
                pp[i].dim[0].org =
                pp[i].dim[1].org = 0.0;
                if (inpres.xr <= inpres.yr*inpaspect) {
                        pp[i].dim[0].siz = inpaspect *
                                                (double)inpres.yr/inpres.xr;
                        pp[i].dim[1].siz = 1.0;
                } else {
                        pp[i].dim[0].siz = 1.0;
                        pp[i].dim[1].siz = (double)inpres.xr/inpres.yr /
                                                inpaspect;
                }
                pp[i].dim[0].p = pp[i].dim[1].p = NULL;
                pp[i].arr = (DATATYPE *)
                                malloc(inpres.xr*inpres.yr*sizeof(DATATYPE));
                if (pp[i].arr == NULL)
                        goto memerr;
        }
                                                        /* load picture */
        sl = scanlen(&inpres);
        ns = numscans(&inpres);
        if ((scanin = (COLOR *)malloc(sl*sizeof(COLOR))) == NULL)
                goto memerr;
        for (y = 0; y < ns; y++) {
                if (freadscan(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);
                        pp[0].arr[i] = colval(scanin[x],RED);
                        pp[1].arr[i] = colval(scanin[x],GRN);
                        pp[2].arr[i] = colval(scanin[x],BLU);
                }
        }
        free((char *)scanin);
        fclose(fp);
        i = hash(pname);
        pp[0].next =
        pp[1].next =
        pp[2].next = ptab[i];
        return(ptab[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;
                                free((char *)dp->arr);
                                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;
        }
}


freepict(pname)                 /* free memory associated with pname */
char  *pname;
{
        DATARRAY  head;
        int  hval, nents;
        register DATARRAY  *pp, *ppl;

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