src/rt/data.c

#ifndef lint
static const char       RCSid[] = "$Id: data.c,v 2.28 2004/10/23 18:55:52 schorsch Exp $";
#endif
/*
 *  data.c - routines dealing with interpolated data.
 */

#include "copyright.h"

#include  <time.h>

#include  "platform.h"
#include  "paths.h"
#include  "standard.h"
#include  "color.h"
#include  "resolu.h"
#include  "view.h"
#include  "data.h"

                                /* picture memory usage before warning */
#ifndef PSIZWARN
#ifdef SMLMEM
#define PSIZWARN        3000000
#else
#define PSIZWARN        10000000
#endif
#endif

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

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


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

static gethfunc headaspect;


extern DATARRAY *
getdata(                                /* 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, getrlibpath(), 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', (char *)&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, (char *)&dp->dim[i].org) <= 0)
                        goto scanerr;
                if (fgetval(fp, DATATY, (char *)&dp->dim[i].siz) <= 0)
                        goto scanerr;
                if (fgetval(fp, 'i', (char *)&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,
                                                (char *)&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, (char *)&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);
        return NULL; /* pro forma return */
}


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

        if (isaspect(s))
                *(double*)iap *= aspectval(s);
        else if (formatval(fmt, s) && !globmatch(PICFMT, fmt))
                *(double*)iap = 0.0;
        return(0);
}


extern DATARRAY *
getpict(                                /* get picture pname */
        char  *pname
)
{
        double  inpaspect;
        char  *pfname;
        FILE  *fp;
        COLR  *scanin;
        int  sl, ns;
        RESOLU  inpres;
        RREAL  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, getrlibpath(), 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);
        }
        SET_FILE_BINARY(fp);
                                                /* 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 %.1f MB of memory",
                                pname, i*(1.0/(1024*1024)));
                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((void *)scanin);
        fclose(fp);
        i = hash(pname);
        pp[0].next = dtab[i];           /* link into picture list */
        pp[1] = pp[0];
        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);
        return NULL; /* pro forma return */
}


extern void
freedata(                       /* release data array reference */
        DATARRAY  *dta
)
{
        DATARRAY  head;
        int  hval, nents;
        register DATARRAY  *dpl, *dp;
        register int  i;

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


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