src/rt/func.c

#ifndef lint
static const char       RCSid[] = "$Id: func.c,v 2.22 2004/03/30 16:13:01 schorsch Exp $";
#endif
/*
 *  func.c - interface to calcomp functions.
 */

#include "copyright.h"

#include <math.h>

#include  "ray.h"
#include  "paths.h"
#include  "otypes.h"
#include  "func.h"


#define  INITFILE       "rayinit.cal"
#define  CALSUF         ".cal"
#define  LCALSUF        4
char  REFVNAME[] = "`FILE_REFCNT";

XF  unitxf = {                  /* identity transform */
        {{1.0, 0.0, 0.0, 0.0},
        {0.0, 1.0, 0.0, 0.0},
        {0.0, 0.0, 1.0, 0.0},
        {0.0, 0.0, 0.0, 1.0}},
        1.0
};

XF  funcxf;                     /* current transformation */
static OBJREC  *fobj = NULL;    /* current function object */
static RAY  *fray = NULL;       /* current function ray */

static double  l_erf(char *), l_erfc(char *), l_arg(char *);


extern MFUNC *
getfunc(        /* get function for this modifier */
        OBJREC  *m,
        int  ff,
        unsigned int  ef,
        int  dofwd
)
{
        static char  initfile[] = INITFILE;
        char  sbuf[MAXSTR];
        register char  **arg;
        register MFUNC  *f;
        int  ne, na;
        register int  i;
                                        /* check to see if done already */
        if ((f = (MFUNC *)m->os) != NULL)
                return(f);
        fobj = NULL; fray = NULL;
        if (initfile[0]) {              /* initialize on first call */
                esupport |= E_VARIABLE|E_FUNCTION|E_INCHAN|E_RCONST|E_REDEFW;
                esupport &= ~(E_OUTCHAN);
                setcontext("");
                scompile("Dx=$1;Dy=$2;Dz=$3;", NULL, 0);
                scompile("Nx=$4;Ny=$5;Nz=$6;", NULL, 0);
                scompile("Px=$7;Py=$8;Pz=$9;", NULL, 0);
                scompile("T=$10;Ts=$25;Rdot=$11;", NULL, 0);
                scompile("S=$12;Tx=$13;Ty=$14;Tz=$15;", NULL, 0);
                scompile("Ix=$16;Iy=$17;Iz=$18;", NULL, 0);
                scompile("Jx=$19;Jy=$20;Jz=$21;", NULL, 0);
                scompile("Kx=$22;Ky=$23;Kz=$24;", NULL, 0);
                scompile("Lu=$26;Lv=$27;", NULL, 0);
                funset("arg", 1, '=', l_arg);
                funset("erf", 1, ':', l_erf);
                funset("erfc", 1, ':', l_erfc);
                setnoisefuncs();
                setprismfuncs();
                loadfunc(initfile);
                initfile[0] = '\0';
        }
        if ((na = m->oargs.nsargs) <= ff)
                goto toofew;
        arg = m->oargs.sarg;
        if ((f = (MFUNC *)calloc(1, sizeof(MFUNC))) == NULL)
                goto memerr;
        i = strlen(arg[ff]);                    /* set up context */
        if (i == 1 && arg[ff][0] == '.')
                setcontext(f->ctx = "");        /* "." means no file */
        else {
                strcpy(sbuf,arg[ff]);           /* file name is context */
                if (i > LCALSUF && !strcmp(sbuf+i-LCALSUF, CALSUF))
                        sbuf[i-LCALSUF] = '\0'; /* remove suffix */
                setcontext(f->ctx = savestr(sbuf));
                if (!vardefined(REFVNAME)) {    /* file loaded? */
                        loadfunc(arg[ff]);
                        varset(REFVNAME, '=', 1.0);
                } else                          /* reference_count++ */
                        varset(REFVNAME, '=', varvalue(REFVNAME)+1.0);
        }
        curfunc = NULL;                 /* parse expressions */
        sprintf(sbuf, "%s \"%s\"", ofun[m->otype].funame, m->oname);
        for (i=0, ne=0; ef && i < na; i++, ef>>=1)
                if (ef & 1) {                   /* flagged as an expression? */
                        if (ne >= MAXEXPR)
                                objerror(m, INTERNAL, "too many expressions");
                        initstr(arg[i], sbuf, 0);
                        f->ep[ne++] = getE1();
                        if (nextc != EOF)
                                syntax("unexpected character");
                }
        if (ef)
                goto toofew;
        if (i <= ff)                    /* find transform args */
                i = ff+1;
        while (i < na && arg[i][0] != '-')
                i++;
        if (i == na)                    /* no transform */
                f->f = f->b = &unitxf;
        else {                          /* get transform */
                if ((f->b = (XF *)malloc(sizeof(XF))) == NULL)
                        goto memerr;
                if (invxf(f->b, na-i, arg+i) != na-i)
                        objerror(m, USER, "bad transform");
                if (f->b->sca < 0.0)
                        f->b->sca = -f->b->sca;
                if (dofwd) {                    /* do both transforms */
                        if ((f->f = (XF *)malloc(sizeof(XF))) == NULL)
                                goto memerr;
                        xf(f->f, na-i, arg+i);
                        if (f->f->sca < 0.0)
                                f->f->sca = -f->f->sca;
                }
        }
        m->os = (char *)f;
        return(f);
toofew:
        objerror(m, USER, "too few string arguments");
memerr:
        error(SYSTEM, "out of memory in getfunc");
        return NULL; /* pro forma return */
}


extern void
freefunc(                       /* free memory associated with modifier */
        OBJREC  *m
)
{
        register MFUNC  *f;
        register int  i;

        if ((f = (MFUNC *)m->os) == NULL)
                return;
        for (i = 0; f->ep[i] != NULL; i++)
                epfree(f->ep[i]);
        if (f->ctx[0]) {                        /* done with definitions */
                setcontext(f->ctx);
                i = varvalue(REFVNAME)-.5;      /* reference_count-- */
                if (i > 0)
                        varset(REFVNAME, '=', (double)i);
                else
                        dcleanup(2);            /* remove definitions */
                freestr(f->ctx);
        }
        if (f->b != &unitxf)
                free((void *)f->b);
        if (f->f != NULL && f->f != &unitxf)
                free((void *)f->f);
        free((void *)f);
        m->os = NULL;
}


extern int
setfunc(                        /* set channels for function call */
        OBJREC  *m,
        register RAY  *r
)
{
        static RNUMBER  lastrno = ~0;
        register MFUNC  *f;
                                        /* get function */
        if ((f = (MFUNC *)m->os) == NULL)
                objerror(m, CONSISTENCY, "setfunc called before getfunc");
                                        /* set evaluator context */
        setcontext(f->ctx);
                                        /* check to see if matrix set */
        if (m == fobj && r->rno == lastrno)
                return(0);
        fobj = m;
        fray = r;
        if (r->rox != NULL)
                if (f->b != &unitxf) {
                        funcxf.sca = r->rox->b.sca * f->b->sca;
                        multmat4(funcxf.xfm, r->rox->b.xfm, f->b->xfm);
                } else
                        funcxf = r->rox->b;
        else
                funcxf = *(f->b);
        lastrno = r->rno;
        eclock++;               /* notify expression evaluator */
        return(1);
}


extern void
loadfunc(                       /* load definition file */
        char  *fname
)
{
        char  *ffname;

        if ((ffname = getpath(fname, getrlibpath(), R_OK)) == NULL) {
                sprintf(errmsg, "cannot find function file \"%s\"", fname);
                error(USER, errmsg);
        }
        fcompile(ffname);
}


static double
l_arg(char *nm)                 /* return nth real argument */
{
        register int  n;

        if (fobj == NULL)
                syntax("arg(n) used in constant expression");

        n = argument(1) + .5;           /* round to integer */

        if (n < 1)
                return(fobj->oargs.nfargs);

        if (n > fobj->oargs.nfargs) {
                sprintf(errmsg, "missing real argument %d", n);
                objerror(fobj, USER, errmsg);
        }
        return(fobj->oargs.farg[n-1]);
}


static double
l_erf(char *nm)                 /* error function */
{
        return(erf(argument(1)));
}


static double
l_erfc(char *nm)                /* cumulative error function */
{
        return(erfc(argument(1)));
}


extern double
chanvalue(                      /* return channel n to calcomp */
        register int  n
)
{
        if (fray == NULL)
                syntax("ray parameter used in constant expression");

        if (--n < 0)
                goto badchan;

        if (n <= 2)                     /* ray direction */

                return( (       fray->rdir[0]*funcxf.xfm[0][n] +
                                fray->rdir[1]*funcxf.xfm[1][n] +
                                fray->rdir[2]*funcxf.xfm[2][n]  )
                         / funcxf.sca );

        if (n <= 5)                     /* surface normal */

                return( (       fray->ron[0]*funcxf.xfm[0][n-3] +
                                fray->ron[1]*funcxf.xfm[1][n-3] +
                                fray->ron[2]*funcxf.xfm[2][n-3] )
                         / funcxf.sca );

        if (n <= 8)                     /* intersection */

                return( fray->rop[0]*funcxf.xfm[0][n-6] +
                                fray->rop[1]*funcxf.xfm[1][n-6] +
                                fray->rop[2]*funcxf.xfm[2][n-6] +
                                             funcxf.xfm[3][n-6] );

        if (n == 9)                     /* total distance */
                return(raydist(fray,PRIMARY) * funcxf.sca);

        if (n == 10)                    /* dot product (range [-1,1]) */
                return( fray->rod <= -1.0 ? -1.0 :
                        fray->rod >= 1.0 ? 1.0 :
                        fray->rod );

        if (n == 11)                    /* scale */
                return(funcxf.sca);

        if (n <= 14)                    /* origin */
                return(funcxf.xfm[3][n-12]);

        if (n <= 17)                    /* i unit vector */
                return(funcxf.xfm[0][n-15] / funcxf.sca);

        if (n <= 20)                    /* j unit vector */
                return(funcxf.xfm[1][n-18] / funcxf.sca);

        if (n <= 23)                    /* k unit vector */
                return(funcxf.xfm[2][n-21] / funcxf.sca);

        if (n == 24)                    /* single ray (shadow) distance */
                return((fray->rot+raydist(fray->parent,SHADOW)) * funcxf.sca);

        if (n <= 26)                    /* local (u,v) coordinates */
                return(fray->uv[n-25]);
badchan:
        error(USER, "illegal channel number");
        return(0.0);
}
Revision:	2.23
Committed:	Sat Dec 12 00:03:42 2009 UTC (14 years, 4 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad4R0
Changes since 2.22:	+2 -2 lines
Log Message:	Created RNUMBER type (#define) in preparation for 64-bit compiling
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.23	static const char RCSid[] = "$Id: func.c,v 2.22 2004/03/30 16:13:01 schorsch Exp $";
3	greg	1.1	#endif
4			/*
5			* func.c - interface to calcomp functions.
6	greg	2.14	*/
7
8	greg	2.15	#include "copyright.h"
9	greg	1.1
10	schorsch	2.22	#include <math.h>
11
12	greg	1.1	#include "ray.h"
13	greg	2.21	#include "paths.h"
14	greg	1.1	#include "otypes.h"
15	greg	2.2	#include "func.h"
16	greg	1.1
17	greg	2.2
18	greg	1.21	#define INITFILE "rayinit.cal"
19	greg	2.2	#define CALSUF ".cal"
20			#define LCALSUF 4
21	greg	2.10	char REFVNAME[] = "`FILE_REFCNT";
22	greg	1.21
23	greg	1.13	XF unitxf = { /* identity transform */
24	greg	2.8	{{1.0, 0.0, 0.0, 0.0},
25			{0.0, 1.0, 0.0, 0.0},
26			{0.0, 0.0, 1.0, 0.0},
27			{0.0, 0.0, 0.0, 1.0}},
28	greg	1.13	1.0
29			};
30
31	greg	1.12	XF funcxf; /* current transformation */
32			static OBJREC fobj = NULL; / current function object */
33			static RAY fray = NULL; / current function ray */
34	greg	1.1
35	greg	2.20	static double l_erf(char ), l_erfc(char ), l_arg(char *);
36	greg	1.1
37	greg	2.2
38	schorsch	2.22	extern MFUNC *
39			getfunc( /* get function for this modifier */
40			OBJREC *m,
41			int ff,
42			unsigned int ef,
43			int dofwd
44			)
45	greg	2.2	{
46			static char initfile[] = INITFILE;
47			char sbuf[MAXSTR];
48			register char **arg;
49			register MFUNC *f;
50			int ne, na;
51			register int i;
52			/* check to see if done already */
53			if ((f = (MFUNC *)m->os) != NULL)
54			return(f);
55	greg	2.7	fobj = NULL; fray = NULL;
56	greg	2.2	if (initfile[0]) { /* initialize on first call */
57	greg	2.14	esupport \|= E_VARIABLE\|E_FUNCTION\|E_INCHAN\|E_RCONST\|E_REDEFW;
58			esupport &= ~(E_OUTCHAN);
59	greg	2.2	setcontext("");
60			scompile("Dx=$1;Dy=$2;Dz=$3;", NULL, 0);
61			scompile("Nx=$4;Ny=$5;Nz=$6;", NULL, 0);
62			scompile("Px=$7;Py=$8;Pz=$9;", NULL, 0);
63	greg	2.6	scompile("T=$10;Ts=$25;Rdot=$11;", NULL, 0);
64	greg	2.2	scompile("S=$12;Tx=$13;Ty=$14;Tz=$15;", NULL, 0);
65			scompile("Ix=$16;Iy=$17;Iz=$18;", NULL, 0);
66			scompile("Jx=$19;Jy=$20;Jz=$21;", NULL, 0);
67			scompile("Kx=$22;Ky=$23;Kz=$24;", NULL, 0);
68	greg	2.17	scompile("Lu=$26;Lv=$27;", NULL, 0);
69	greg	2.2	funset("arg", 1, '=', l_arg);
70			funset("erf", 1, ':', l_erf);
71			funset("erfc", 1, ':', l_erfc);
72			setnoisefuncs();
73	greg	2.11	setprismfuncs();
74	greg	2.2	loadfunc(initfile);
75			initfile[0] = '\0';
76			}
77			if ((na = m->oargs.nsargs) <= ff)
78			goto toofew;
79			arg = m->oargs.sarg;
80			if ((f = (MFUNC *)calloc(1, sizeof(MFUNC))) == NULL)
81			goto memerr;
82			i = strlen(arg[ff]); /* set up context */
83			if (i == 1 && arg[ff][0] == '.')
84			setcontext(f->ctx = ""); /* "." means no file */
85			else {
86	greg	2.14	strcpy(sbuf,arg[ff]); /* file name is context */
87	greg	2.2	if (i > LCALSUF && !strcmp(sbuf+i-LCALSUF, CALSUF))
88			sbuf[i-LCALSUF] = '\0'; /* remove suffix */
89			setcontext(f->ctx = savestr(sbuf));
90			if (!vardefined(REFVNAME)) { /* file loaded? */
91	greg	2.5	loadfunc(arg[ff]);
92	greg	2.2	varset(REFVNAME, '=', 1.0);
93			} else /* reference_count++ */
94			varset(REFVNAME, '=', varvalue(REFVNAME)+1.0);
95			}
96			curfunc = NULL; /* parse expressions */
97			sprintf(sbuf, "%s \"%s\"", ofun[m->otype].funame, m->oname);
98			for (i=0, ne=0; ef && i < na; i++, ef>>=1)
99			if (ef & 1) { /* flagged as an expression? */
100			if (ne >= MAXEXPR)
101			objerror(m, INTERNAL, "too many expressions");
102			initstr(arg[i], sbuf, 0);
103			f->ep[ne++] = getE1();
104			if (nextc != EOF)
105			syntax("unexpected character");
106			}
107			if (ef)
108			goto toofew;
109			if (i <= ff) /* find transform args */
110			i = ff+1;
111			while (i < na && arg[i][0] != '-')
112			i++;
113			if (i == na) /* no transform */
114			f->f = f->b = &unitxf;
115			else { /* get transform */
116			if ((f->b = (XF *)malloc(sizeof(XF))) == NULL)
117	greg	2.3	goto memerr;
118	greg	2.2	if (invxf(f->b, na-i, arg+i) != na-i)
119			objerror(m, USER, "bad transform");
120			if (f->b->sca < 0.0)
121			f->b->sca = -f->b->sca;
122			if (dofwd) { /* do both transforms */
123			if ((f->f = (XF *)malloc(sizeof(XF))) == NULL)
124	greg	2.3	goto memerr;
125	greg	2.2	xf(f->f, na-i, arg+i);
126			if (f->f->sca < 0.0)
127			f->f->sca = -f->f->sca;
128			}
129			}
130			m->os = (char *)f;
131			return(f);
132			toofew:
133	greg	2.4	objerror(m, USER, "too few string arguments");
134	greg	2.2	memerr:
135			error(SYSTEM, "out of memory in getfunc");
136	schorsch	2.22	return NULL; /* pro forma return */
137	greg	2.2	}
138
139
140	schorsch	2.22	extern void
141			freefunc( /* free memory associated with modifier */
142			OBJREC *m
143			)
144	greg	2.2	{
145			register MFUNC *f;
146			register int i;
147
148			if ((f = (MFUNC *)m->os) == NULL)
149			return;
150			for (i = 0; f->ep[i] != NULL; i++)
151			epfree(f->ep[i]);
152			if (f->ctx[0]) { /* done with definitions */
153			setcontext(f->ctx);
154			i = varvalue(REFVNAME)-.5; /* reference_count-- */
155			if (i > 0)
156			varset(REFVNAME, '=', (double)i);
157			else
158			dcleanup(2); /* remove definitions */
159			freestr(f->ctx);
160			}
161			if (f->b != &unitxf)
162	greg	2.14	free((void *)f->b);
163	greg	2.2	if (f->f != NULL && f->f != &unitxf)
164	greg	2.14	free((void *)f->f);
165			free((void *)f);
166	greg	2.2	m->os = NULL;
167			}
168
169
170	schorsch	2.22	extern int
171			setfunc( /* set channels for function call */
172			OBJREC *m,
173			register RAY *r
174			)
175	greg	1.1	{
176	greg	2.23	static RNUMBER lastrno = ~0;
177	greg	2.2	register MFUNC *f;
178			/* get function */
179			if ((f = (MFUNC *)m->os) == NULL)
180	greg	2.14	objerror(m, CONSISTENCY, "setfunc called before getfunc");
181	greg	2.2	/* set evaluator context */
182			setcontext(f->ctx);
183			/* check to see if matrix set */
184	greg	1.16	if (m == fobj && r->rno == lastrno)
185	greg	1.17	return(0);
186	greg	1.1	fobj = m;
187			fray = r;
188	greg	1.13	if (r->rox != NULL)
189	greg	2.2	if (f->b != &unitxf) {
190			funcxf.sca = r->rox->b.sca * f->b->sca;
191			multmat4(funcxf.xfm, r->rox->b.xfm, f->b->xfm);
192	greg	1.13	} else
193	schorsch	2.19	funcxf = r->rox->b;
194	greg	1.13	else
195	schorsch	2.19	funcxf = *(f->b);
196	greg	1.16	lastrno = r->rno;
197	greg	1.1	eclock++; /* notify expression evaluator */
198	greg	1.17	return(1);
199	greg	1.1	}
200
201
202	schorsch	2.22	extern void
203			loadfunc( /* load definition file */
204			char *fname
205			)
206	greg	1.1	{
207			char *ffname;
208
209	greg	2.18	if ((ffname = getpath(fname, getrlibpath(), R_OK)) == NULL) {
210	greg	1.1	sprintf(errmsg, "cannot find function file \"%s\"", fname);
211			error(USER, errmsg);
212			}
213			fcompile(ffname);
214			}
215
216
217	greg	2.2	static double
218	greg	2.20	l_arg(char nm) / return nth real argument */
219	greg	1.1	{
220			register int n;
221
222	greg	2.7	if (fobj == NULL)
223			syntax("arg(n) used in constant expression");
224
225	greg	1.1	n = argument(1) + .5; /* round to integer */
226
227			if (n < 1)
228			return(fobj->oargs.nfargs);
229
230			if (n > fobj->oargs.nfargs) {
231			sprintf(errmsg, "missing real argument %d", n);
232			objerror(fobj, USER, errmsg);
233			}
234			return(fobj->oargs.farg[n-1]);
235	greg	1.19	}
236
237
238	greg	2.2	static double
239	greg	2.20	l_erf(char nm) / error function */
240	greg	1.19	{
241			return(erf(argument(1)));
242			}
243
244
245	greg	2.2	static double
246	greg	2.20	l_erfc(char nm) / cumulative error function */
247	greg	1.19	{
248			return(erfc(argument(1)));
249	greg	1.1	}
250
251
252	schorsch	2.22	extern double
253			chanvalue( /* return channel n to calcomp */
254			register int n
255			)
256	greg	1.1	{
257	greg	2.7	if (fray == NULL)
258			syntax("ray parameter used in constant expression");
259	greg	1.1
260	greg	1.10	if (--n < 0)
261			goto badchan;
262	greg	1.1
263	greg	2.17	if (n <= 2) /* ray direction */
264	greg	1.1
265	greg	1.12	return( ( fray->rdir[0]*funcxf.xfm[0][n] +
266			fray->rdir[1]*funcxf.xfm[1][n] +
267			fray->rdir[2]*funcxf.xfm[2][n] )
268			/ funcxf.sca );
269	greg	1.1
270	greg	2.17	if (n <= 5) /* surface normal */
271	greg	1.8
272	greg	1.12	return( ( fray->ron[0]*funcxf.xfm[0][n-3] +
273			fray->ron[1]*funcxf.xfm[1][n-3] +
274			fray->ron[2]*funcxf.xfm[2][n-3] )
275			/ funcxf.sca );
276	greg	1.8
277	greg	2.17	if (n <= 8) /* intersection */
278	greg	1.8
279	greg	1.12	return( fray->rop[0]*funcxf.xfm[0][n-6] +
280			fray->rop[1]*funcxf.xfm[1][n-6] +
281			fray->rop[2]*funcxf.xfm[2][n-6] +
282			funcxf.xfm[3][n-6] );
283	greg	1.8
284	greg	2.13	if (n == 9) /* total distance */
285			return(raydist(fray,PRIMARY) * funcxf.sca);
286	greg	1.10
287	greg	1.18	if (n == 10) /* dot product (range [-1,1]) */
288			return( fray->rod <= -1.0 ? -1.0 :
289			fray->rod >= 1.0 ? 1.0 :
290			fray->rod );
291	greg	1.10
292	greg	1.8	if (n == 11) /* scale */
293	greg	1.12	return(funcxf.sca);
294	greg	1.8
295	greg	2.17	if (n <= 14) /* origin */
296	greg	1.12	return(funcxf.xfm[3][n-12]);
297	greg	1.8
298	greg	2.17	if (n <= 17) /* i unit vector */
299	greg	1.12	return(funcxf.xfm[0][n-15] / funcxf.sca);
300	greg	1.8
301	greg	2.17	if (n <= 20) /* j unit vector */
302			return(funcxf.xfm[1][n-18] / funcxf.sca);
303	greg	1.8
304	greg	2.17	if (n <= 23) /* k unit vector */
305	greg	1.12	return(funcxf.xfm[2][n-21] / funcxf.sca);
306	greg	2.6
307	greg	2.17	if (n == 24) /* single ray (shadow) distance */
308	greg	2.13	return((fray->rot+raydist(fray->parent,SHADOW)) * funcxf.sca);
309	greg	2.16
310	greg	2.17	if (n <= 26) /* local (u,v) coordinates */
311			return(fray->uv[n-25]);
312	greg	1.10	badchan:
313			error(USER, "illegal channel number");
314	greg	2.16	return(0.0);
315	greg	1.1	}