src/rt/data.c

/* Copyright (c) 1986 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  "data.h"


extern char  *libpath;                  /* library search path */

static DATARRAY  *dlist = NULL;         /* data array list */

static DATARRAY  *plist = NULL;         /* 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 = dlist; 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 directorys 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:
         *
         *              @ ni p0i p1i .. pni
         */

        if ((dfname = getpath(dname, libpath, 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 (fscanf(fp, "%d", &dp->nd) != 1)
                goto scanerr;
        if (dp->nd <= 0 || dp->nd > MAXDIM) {
                sprintf(errmsg, "bad number of dimensions for \"%s\"", dname);
                error(USER, errmsg);
        }
        asize = 1;
        for (i = 0; i < dp->nd; i++) {
                if (fscanf(fp, "%lf %lf %d",
                                &dp->dim[i].org, &dp->dim[i].siz,
                                &dp->dim[i].ne) == 3) {
                        dp->dim[i].siz -= dp->dim[i].org;
                        dp->dim[i].p = NULL;
                } else if (fscanf(fp, "@ %d", &dp->dim[i].ne) == 1) {
                        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 (fscanf(fp, "%lf", &dp->dim[i].p[j]) != 1)
                                        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
                        goto scanerr;
                if (dp->dim[i].siz == 0.0 || dp->dim[i].ne < 2)
                        goto scanerr;
                asize *= dp->dim[i].ne;
        }
        if ((dp->arr = (DATATYPE *)malloc(asize*sizeof(DATATYPE))) == NULL)
                goto memerr;
        
        for (i = 0; i < asize; i++)
                if (fscanf(fp, DSCANF, &dp->arr[i]) != 1)
                        goto scanerr;
        
        fclose(fp);
        dp->next = dlist;
        return(dlist = 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);
}


DATARRAY *
getpict(pname)                          /* get picture pname */
char  *pname;
{
        extern char  *libpath;
        char  *pfname;
        FILE  *fp;
        COLOR  *scanin;
        int  width, height;
        int  x, y;
        register int  i;
        register DATARRAY  *pp;
                                                /* look for array in list */
        for (pp = plist; pp != NULL; pp = pp->next)
                if (!strcmp(pname, pp->name))
                        return(pp);             /* found! */

        if ((pfname = getpath(pname, libpath, 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);
        }
                                                /* get dimensions */
        getheader(fp, NULL);
        if (fgetresolu(&width, &height, fp) != (YMAJOR|YDECR))
                goto readerr;
        for (i = 0; i < 3; i++) {
                pp[i].nd = 2;
                pp[i].dim[0].ne = width;
                pp[i].dim[1].ne = height;
                pp[i].dim[0].org =
                pp[i].dim[1].org = 0.0;
                if (width <= height) {
                        pp[i].dim[0].siz = 1.0;
                        pp[i].dim[1].siz = (double)height/width;
                } else {
                        pp[i].dim[0].siz = (double)width/height;
                        pp[i].dim[1].siz = 1.0;
                }
                pp[i].arr = (DATATYPE *)malloc(width*height*sizeof(DATATYPE));
                if (pp[i].arr == NULL)
                        goto memerr;
        }
                                                        /* load picture */
        if ((scanin = (COLOR *)malloc(width*sizeof(COLOR))) == NULL)
                goto memerr;
        for (y = height-1; y >= 0; y--) {
                if (freadscan(scanin, width, fp) < 0)
                        goto readerr;
                for (x = 0; x < width; x++)
                        for (i = 0; i < 3; i++)
                                pp[i].arr[x*height+y] = colval(scanin[x],i);
        }
        free((char *)scanin);
        fclose(fp);
        pp[0].next =
        pp[1].next =
        pp[2].next = plist;
        return(plist = 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;
{
        register DATARRAY  *dp, *dpl;
        register int  i;

        for (dpl = NULL, dp = dlist; dp != NULL; dpl = dp, dp = dp->next)
                if (!strcmp(dname, dp->name)) {
                        if (dpl == NULL)
                                dlist = dp->next;
                        else
                                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);
                        return;
                }
}


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

        for (ppl = NULL, pp = plist; pp != NULL; ppl = pp, pp = pp->next)
                if (!strcmp(pname, pp->name)) {
                        if (ppl == NULL)
                                plist = pp->next;
                        else
                                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);
                        return;
                }
}


double
datavalue(dp, pt)               /* interpolate data value at a point */
register DATARRAY  *dp;
double  *pt;
{
        DATARRAY  sd;
        int  asize;
        register int  i;
        double  x, y, 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 = x * (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)
                        for (i = 0; i < dp->dim[0].ne; i++)
                                if (pt[0] < dp->dim[0].p[i])
                                        break;
                else
                        for (i = 0; i < dp->dim[0].ne; i++)
                                if (pt[0] >= dp->dim[0].p[i])
                                        break;
                if (i <= 0)
                        i = 0;
                else if (i >= dp->dim[0].ne)
                        i = dp->dim[0].ne - 2;
                else
                        i--;
                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)
                y = (2*y1-y0)/(x-i-1);
        else if (x < i-1)
                y = (2*y0-y1)/(i-x);
        else
                y = y0*((i+1)-x) + y1*(x-i);

        return(y);
}
Revision:	1.4
Committed:	Tue Apr 3 17:24:35 1990 UTC (34 years, 1 month ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.3:	+61 -17 lines
Log Message:	Added capability for unevenly spaced points
#	User	Rev	Content
1	greg	1.1	/* Copyright (c) 1986 Regents of the University of California */
2
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			#include "data.h"
18
19
20			extern char libpath; / library search path */
21
22			static DATARRAY dlist = NULL; / data array list */
23
24			static DATARRAY plist = NULL; / picture list */
25
26
27			DATARRAY *
28			getdata(dname) /* get data array dname */
29			char *dname;
30			{
31			char *dfname;
32			FILE *fp;
33			int asize;
34	greg	1.4	register int i, j;
35	greg	1.1	register DATARRAY *dp;
36			/* look for array in list */
37			for (dp = dlist; dp != NULL; dp = dp->next)
38			if (!strcmp(dname, dp->name))
39			return(dp); /* found! */
40
41			/*
42			* If we haven't loaded the data already, we will look
43			* for it in the directorys specified by the library path.
44			*
45			* The file has the following format:
46			*
47	greg	1.4	* N
48	greg	1.1	* beg0 end0 n0
49			* beg1 end1 n1
50			* . . .
51	greg	1.4	* begN endN nN
52	greg	1.1	* data, later dimensions changing faster
53			* . . .
54			*
55	greg	1.4	* For irregularly spaced points, the following can be
56			* substituted for begi endi ni:
57			*
58			* @ ni p0i p1i .. pni
59	greg	1.1	*/
60
61	greg	1.2	if ((dfname = getpath(dname, libpath, R_OK)) == NULL) {
62	greg	1.1	sprintf(errmsg, "cannot find data file \"%s\"", dname);
63			error(USER, errmsg);
64			}
65			if ((dp = (DATARRAY *)malloc(sizeof(DATARRAY))) == NULL)
66			goto memerr;
67
68			dp->name = savestr(dname);
69
70			if ((fp = fopen(dfname, "r")) == NULL) {
71			sprintf(errmsg, "cannot open data file \"%s\"", dfname);
72			error(SYSTEM, errmsg);
73			}
74			/* get dimensions */
75			if (fscanf(fp, "%d", &dp->nd) != 1)
76			goto scanerr;
77			if (dp->nd <= 0 \|\| dp->nd > MAXDIM) {
78			sprintf(errmsg, "bad number of dimensions for \"%s\"", dname);
79			error(USER, errmsg);
80			}
81			asize = 1;
82			for (i = 0; i < dp->nd; i++) {
83			if (fscanf(fp, "%lf %lf %d",
84			&dp->dim[i].org, &dp->dim[i].siz,
85	greg	1.4	&dp->dim[i].ne) == 3) {
86			dp->dim[i].siz -= dp->dim[i].org;
87			dp->dim[i].p = NULL;
88			} else if (fscanf(fp, "@ %d", &dp->dim[i].ne) == 1) {
89			dp->dim[i].p = (double )malloc(dp->dim[i].nesizeof(double));
90			if (dp->dim[i].p == NULL)
91			goto memerr;
92			for (j = 0; j < dp->dim[i].ne; j++)
93			if (fscanf(fp, "%lf", &dp->dim[i].p[j]) != 1)
94			goto scanerr;
95			for (j = 1; j < dp->dim[i].ne-1; j++)
96			if ((dp->dim[i].p[j-1] < dp->dim[i].p[j]) !=
97			(dp->dim[i].p[j] < dp->dim[i].p[j+1]))
98			goto scanerr;
99			dp->dim[i].org = dp->dim[i].p[0];
100			dp->dim[i].siz = dp->dim[i].p[dp->dim[i].ne-1]
101			- dp->dim[i].p[0];
102			} else
103	greg	1.1	goto scanerr;
104	greg	1.4	if (dp->dim[i].siz == 0.0 \|\| dp->dim[i].ne < 2)
105			goto scanerr;
106	greg	1.1	asize *= dp->dim[i].ne;
107			}
108			if ((dp->arr = (DATATYPE )malloc(asizesizeof(DATATYPE))) == NULL)
109			goto memerr;
110
111			for (i = 0; i < asize; i++)
112			if (fscanf(fp, DSCANF, &dp->arr[i]) != 1)
113			goto scanerr;
114
115			fclose(fp);
116			dp->next = dlist;
117			return(dlist = dp);
118
119			memerr:
120			error(SYSTEM, "out of memory in getdata");
121			scanerr:
122			sprintf(errmsg, "%s in data file \"%s\"",
123			feof(fp) ? "unexpected EOF" : "bad format", dfname);
124			error(USER, errmsg);
125			}
126
127
128			DATARRAY *
129			getpict(pname) /* get picture pname */
130			char *pname;
131			{
132			extern char *libpath;
133			char *pfname;
134			FILE *fp;
135			COLOR *scanin;
136			int width, height;
137			int x, y;
138			register int i;
139			register DATARRAY *pp;
140			/* look for array in list */
141			for (pp = plist; pp != NULL; pp = pp->next)
142			if (!strcmp(pname, pp->name))
143			return(pp); /* found! */
144
145	greg	1.2	if ((pfname = getpath(pname, libpath, R_OK)) == NULL) {
146	greg	1.1	sprintf(errmsg, "cannot find picture file \"%s\"", pname);
147			error(USER, errmsg);
148			}
149			if ((pp = (DATARRAY *)calloc(3, sizeof(DATARRAY))) == NULL)
150			goto memerr;
151
152			pp[0].name =
153			pp[1].name =
154			pp[2].name = savestr(pname);
155
156			if ((fp = fopen(pfname, "r")) == NULL) {
157			sprintf(errmsg, "cannot open picture file \"%s\"", pfname);
158			error(SYSTEM, errmsg);
159			}
160			/* get dimensions */
161			getheader(fp, NULL);
162	greg	1.3	if (fgetresolu(&width, &height, fp) != (YMAJOR\|YDECR))
163	greg	1.1	goto readerr;
164			for (i = 0; i < 3; i++) {
165			pp[i].nd = 2;
166			pp[i].dim[0].ne = width;
167			pp[i].dim[1].ne = height;
168			pp[i].dim[0].org =
169			pp[i].dim[1].org = 0.0;
170			if (width <= height) {
171			pp[i].dim[0].siz = 1.0;
172			pp[i].dim[1].siz = (double)height/width;
173			} else {
174			pp[i].dim[0].siz = (double)width/height;
175			pp[i].dim[1].siz = 1.0;
176			}
177			pp[i].arr = (DATATYPE )malloc(widthheight*sizeof(DATATYPE));
178			if (pp[i].arr == NULL)
179			goto memerr;
180			}
181			/* load picture */
182			if ((scanin = (COLOR )malloc(widthsizeof(COLOR))) == NULL)
183			goto memerr;
184			for (y = height-1; y >= 0; y--) {
185			if (freadscan(scanin, width, fp) < 0)
186			goto readerr;
187			for (x = 0; x < width; x++)
188			for (i = 0; i < 3; i++)
189			pp[i].arr[x*height+y] = colval(scanin[x],i);
190			}
191			free((char *)scanin);
192			fclose(fp);
193			pp[0].next =
194			pp[1].next =
195			pp[2].next = plist;
196			return(plist = pp);
197
198			memerr:
199			error(SYSTEM, "out of memory in getpict");
200			readerr:
201			sprintf(errmsg, "bad picture file \"%s\"", pfname);
202			error(USER, errmsg);
203			}
204
205
206			freedata(dname) /* free memory associated with dname */
207			char *dname;
208			{
209			register DATARRAY dp, dpl;
210	greg	1.4	register int i;
211	greg	1.1
212			for (dpl = NULL, dp = dlist; dp != NULL; dpl = dp, dp = dp->next)
213			if (!strcmp(dname, dp->name)) {
214			if (dpl == NULL)
215			dlist = dp->next;
216			else
217			dpl->next = dp->next;
218			free((char *)dp->arr);
219	greg	1.4	for (i = 0; i < dp->nd; i++)
220			if (dp->dim[i].p != NULL)
221			free((char *)dp->dim[i].p);
222	greg	1.1	freestr(dp->name);
223			free((char *)dp);
224			return;
225			}
226			}
227
228
229			freepict(pname) /* free memory associated with pname */
230			char *pname;
231			{
232			register DATARRAY pp, ppl;
233
234			for (ppl = NULL, pp = plist; pp != NULL; ppl = pp, pp = pp->next)
235			if (!strcmp(pname, pp->name)) {
236			if (ppl == NULL)
237			plist = pp->next;
238			else
239			ppl->next = pp->next;
240			free((char *)pp[0].arr);
241			free((char *)pp[1].arr);
242			free((char *)pp[2].arr);
243			freestr(pp[0].name);
244			free((char *)pp);
245			return;
246			}
247			}
248
249
250			double
251			datavalue(dp, pt) /* interpolate data value at a point */
252			register DATARRAY *dp;
253			double *pt;
254			{
255			DATARRAY sd;
256			int asize;
257			register int i;
258			double x, y, y0, y1;
259	greg	1.4	/* set up dimensions for recursion */
260	greg	1.1	sd.nd = dp->nd - 1;
261			asize = 1;
262			for (i = 0; i < sd.nd; i++) {
263			sd.dim[i].org = dp->dim[i+1].org;
264			sd.dim[i].siz = dp->dim[i+1].siz;
265	greg	1.4	sd.dim[i].p = dp->dim[i+1].p;
266	greg	1.1	asize *= sd.dim[i].ne = dp->dim[i+1].ne;
267			}
268	greg	1.4	/* get independent variable */
269			if (dp->dim[0].p == NULL) { /* evenly spaced points */
270			x = (pt[0] - dp->dim[0].org)/dp->dim[0].siz;
271			x = x * (dp->dim[0].ne - 1);
272			i = x;
273			if (i < 0)
274			i = 0;
275			else if (i > dp->dim[0].ne - 2)
276			i = dp->dim[0].ne - 2;
277			} else { /* unevenly spaced points */
278			if (dp->dim[0].siz > 0.0)
279			for (i = 0; i < dp->dim[0].ne; i++)
280			if (pt[0] < dp->dim[0].p[i])
281			break;
282			else
283			for (i = 0; i < dp->dim[0].ne; i++)
284			if (pt[0] >= dp->dim[0].p[i])
285			break;
286			if (i <= 0)
287			i = 0;
288			else if (i >= dp->dim[0].ne)
289			i = dp->dim[0].ne - 2;
290			else
291			i--;
292			x = i + (pt[0] - dp->dim[0].p[i]) /
293			(dp->dim[0].p[i+1] - dp->dim[0].p[i]);
294			}
295			/* get dependent variable */
296	greg	1.1	if (dp->nd == 1) {
297			y0 = dp->arr[i];
298			y1 = dp->arr[i+1];
299			} else {
300			sd.arr = &dp->arr[i*asize];
301			y0 = datavalue(&sd, pt+1);
302			sd.arr = &dp->arr[(i+1)*asize];
303			y1 = datavalue(&sd, pt+1);
304			}
305			/*
306			* Extrapolate as far as one division, then
307			* taper off harmonically to zero.
308			*/
309			if (x > i+2)
310			y = (2*y1-y0)/(x-i-1);
311			else if (x < i-1)
312			y = (2*y0-y1)/(i-x);
313			else
314			y = y0((i+1)-x) + y1(x-i);
315
316			return(y);
317			}