src/rt/func.c

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

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

/*
 *  func.c - interface to calcomp functions.
 *
 *     4/7/86
 */

#include  "ray.h"

#include  "otypes.h"


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 */


setmap(m, r, bx)                /* set channels for function call */
OBJREC  *m;
register RAY  *r;
XF  *bx;
{
        extern double  l_arg(), l_erf(), l_erfc();
        extern long  eclock;
        static char  *initfile = "rayinit.cal";
        static long  lastrno = -1;
                                        /* check to see if already set */
        if (m == fobj && r->rno == lastrno)
                return(0);
                                        /* initialize if first call */
        if (initfile != NULL) {
                loadfunc(initfile);
                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;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);
                funset("arg", 1, '=', l_arg);
                funset("erf", 1, ':', l_erf);
                funset("erfc", 1, ':', l_erfc);
                setnoisefuncs();
                initfile = NULL;
        }
        fobj = m;
        fray = r;
        if (r->rox != NULL)
                if (bx != &unitxf) {
                        funcxf.sca = r->rox->b.sca * bx->sca;
                        multmat4(funcxf.xfm, r->rox->b.xfm, bx->xfm);
                } else
                        copystruct(&funcxf, &r->rox->b);
        else
                copystruct(&funcxf, bx);
        lastrno = r->rno;
        eclock++;               /* notify expression evaluator */
        return(1);
}


setfunc(m, r)                           /* simplified interface to setmap */
register OBJREC  *m;
RAY  *r;
{
        register XF  *mxf;

        if ((mxf = (XF *)m->os) == NULL) {
                register int  n;
                register char  **sa;

                for (n = m->oargs.nsargs, sa = m->oargs.sarg;
                                n > 0 && **sa != '-'; n--, sa++)
                        ;
                if (n == 0)
                        mxf = &unitxf;
                else {
                        mxf = (XF *)malloc(sizeof(XF));
                        if (mxf == NULL)
                                goto memerr;
                        if (invxf(mxf, n, sa) != n)
                                objerror(m, USER, "bad transform");
                        if (mxf->sca < 0.0)
                                mxf->sca = -mxf->sca;
                }
                m->os = (char *)mxf;
        }
        return(setmap(m, r, mxf));
memerr:
        error(SYSTEM, "out of memory in setfunc");
}


loadfunc(fname)                 /* load definition file */
char  *fname;
{
        extern char  *libpath;          /* library search path */
        char  *ffname;

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


double
l_arg()                         /* return nth real argument */
{
        extern double  argument();
        register int  n;

        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]);
}


double
l_erf()                         /* error function */
{
        extern double  erf();

        return(erf(argument(1)));
}


double
l_erfc()                        /* cumulative error function */
{
        extern double  erfc();

        return(erfc(argument(1)));
}


double
chanvalue(n)                    /* return channel n to calcomp */
register int  n;
{
        double  sum;
        register RAY  *r;

        if (--n < 0)
                goto badchan;

        if (n < 3)                      /* 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 < 6)                      /* 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 < 9)                      /* 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) {                   /* distance */

                sum = fray->rot;
                for (r = fray->parent; r != NULL; r = r->parent)
                        sum += r->rot;
                return(sum * 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 < 15)                     /* origin */
                return(funcxf.xfm[3][n-12]);

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

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

        if (n < 24)                     /* k unit vector */
                return(funcxf.xfm[2][n-21] / funcxf.sca);
badchan:
        error(USER, "illegal channel number");
}
Revision:	1.19
Committed:	Wed Jun 19 15:08:39 1991 UTC (34 years, 4 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.18:	+21 -1 lines
Log Message:	added erf() and erfc() to function complement
#	User	Rev	Content
1	greg	1.13	/* Copyright (c) 1991 Regents of the University of California */
2	greg	1.1
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6
7			/*
8			* func.c - interface to calcomp functions.
9			*
10			* 4/7/86
11			*/
12
13			#include "ray.h"
14
15			#include "otypes.h"
16
17
18	greg	1.13	XF unitxf = { /* identity transform */
19			1.0, 0.0, 0.0, 0.0,
20			0.0, 1.0, 0.0, 0.0,
21			0.0, 0.0, 1.0, 0.0,
22			0.0, 0.0, 0.0, 1.0,
23			1.0
24			};
25
26	greg	1.12	XF funcxf; /* current transformation */
27			static OBJREC fobj = NULL; / current function object */
28			static RAY fray = NULL; / current function ray */
29	greg	1.1
30
31	greg	1.12	setmap(m, r, bx) /* set channels for function call */
32	greg	1.1	OBJREC *m;
33			register RAY *r;
34	greg	1.12	XF *bx;
35	greg	1.1	{
36	greg	1.19	extern double l_arg(), l_erf(), l_erfc();
37	greg	1.1	extern long eclock;
38			static char *initfile = "rayinit.cal";
39	greg	1.16	static long lastrno = -1;
40			/* check to see if already set */
41			if (m == fobj && r->rno == lastrno)
42	greg	1.17	return(0);
43	greg	1.12	/* initialize if first call */
44	greg	1.1	if (initfile != NULL) {
45			loadfunc(initfile);
46	greg	1.9	scompile("Dx=$1;Dy=$2;Dz=$3;", NULL, 0);
47			scompile("Nx=$4;Ny=$5;Nz=$6;", NULL, 0);
48			scompile("Px=$7;Py=$8;Pz=$9;", NULL, 0);
49			scompile("T=$10;Rdot=$11;", NULL, 0);
50			scompile("S=$12;Tx=$13;Ty=$14;Tz=$15;", NULL, 0);
51			scompile("Ix=$16;Iy=$17;Iz=$18;", NULL, 0);
52			scompile("Jx=$19;Jy=$20;Jz=$21;", NULL, 0);
53			scompile("Kx=$22;Ky=$23;Kz=$24;", NULL, 0);
54	greg	1.14	funset("arg", 1, '=', l_arg);
55	greg	1.19	funset("erf", 1, ':', l_erf);
56			funset("erfc", 1, ':', l_erfc);
57	greg	1.15	setnoisefuncs();
58	greg	1.1	initfile = NULL;
59			}
60			fobj = m;
61			fray = r;
62	greg	1.13	if (r->rox != NULL)
63			if (bx != &unitxf) {
64			funcxf.sca = r->rox->b.sca * bx->sca;
65			multmat4(funcxf.xfm, r->rox->b.xfm, bx->xfm);
66			} else
67			copystruct(&funcxf, &r->rox->b);
68			else
69	greg	1.12	copystruct(&funcxf, bx);
70	greg	1.16	lastrno = r->rno;
71	greg	1.1	eclock++; /* notify expression evaluator */
72	greg	1.17	return(1);
73	greg	1.1	}
74
75
76			setfunc(m, r) /* simplified interface to setmap */
77			register OBJREC *m;
78			RAY *r;
79			{
80	greg	1.2	register XF *mxf;
81	greg	1.1
82	greg	1.2	if ((mxf = (XF *)m->os) == NULL) {
83	greg	1.5	register int n;
84			register char **sa;
85	greg	1.1
86	greg	1.5	for (n = m->oargs.nsargs, sa = m->oargs.sarg;
87			n > 0 && **sa != '-'; n--, sa++)
88			;
89	greg	1.13	if (n == 0)
90			mxf = &unitxf;
91			else {
92			mxf = (XF *)malloc(sizeof(XF));
93			if (mxf == NULL)
94			goto memerr;
95			if (invxf(mxf, n, sa) != n)
96			objerror(m, USER, "bad transform");
97			if (mxf->sca < 0.0)
98			mxf->sca = -mxf->sca;
99			}
100	greg	1.2	m->os = (char *)mxf;
101	greg	1.1	}
102	greg	1.17	return(setmap(m, r, mxf));
103	greg	1.1	memerr:
104			error(SYSTEM, "out of memory in setfunc");
105			}
106
107
108			loadfunc(fname) /* load definition file */
109			char *fname;
110			{
111			extern char libpath; / library search path */
112			char *ffname;
113
114	greg	1.6	if ((ffname = getpath(fname, libpath, R_OK)) == NULL) {
115	greg	1.1	sprintf(errmsg, "cannot find function file \"%s\"", fname);
116			error(USER, errmsg);
117			}
118			fcompile(ffname);
119			}
120
121
122			double
123			l_arg() /* return nth real argument */
124			{
125			extern double argument();
126			register int n;
127
128			n = argument(1) + .5; /* round to integer */
129
130			if (n < 1)
131			return(fobj->oargs.nfargs);
132
133			if (n > fobj->oargs.nfargs) {
134			sprintf(errmsg, "missing real argument %d", n);
135			objerror(fobj, USER, errmsg);
136			}
137			return(fobj->oargs.farg[n-1]);
138	greg	1.19	}
139
140
141			double
142			l_erf() /* error function */
143			{
144			extern double erf();
145
146			return(erf(argument(1)));
147			}
148
149
150			double
151			l_erfc() /* cumulative error function */
152			{
153			extern double erfc();
154
155			return(erfc(argument(1)));
156	greg	1.1	}
157
158
159			double
160			chanvalue(n) /* return channel n to calcomp */
161			register int n;
162			{
163	greg	1.8	double sum;
164	greg	1.1	register RAY *r;
165
166	greg	1.10	if (--n < 0)
167			goto badchan;
168	greg	1.1
169	greg	1.8	if (n < 3) /* ray direction */
170	greg	1.1
171	greg	1.12	return( ( fray->rdir[0]*funcxf.xfm[0][n] +
172			fray->rdir[1]*funcxf.xfm[1][n] +
173			fray->rdir[2]*funcxf.xfm[2][n] )
174			/ funcxf.sca );
175	greg	1.1
176	greg	1.8	if (n < 6) /* surface normal */
177
178	greg	1.12	return( ( fray->ron[0]*funcxf.xfm[0][n-3] +
179			fray->ron[1]*funcxf.xfm[1][n-3] +
180			fray->ron[2]*funcxf.xfm[2][n-3] )
181			/ funcxf.sca );
182	greg	1.8
183			if (n < 9) /* intersection */
184
185	greg	1.12	return( fray->rop[0]*funcxf.xfm[0][n-6] +
186			fray->rop[1]*funcxf.xfm[1][n-6] +
187			fray->rop[2]*funcxf.xfm[2][n-6] +
188			funcxf.xfm[3][n-6] );
189	greg	1.8
190	greg	1.10	if (n == 9) { /* distance */
191
192			sum = fray->rot;
193			for (r = fray->parent; r != NULL; r = r->parent)
194			sum += r->rot;
195	greg	1.12	return(sum * funcxf.sca);
196	greg	1.10
197			}
198	greg	1.18	if (n == 10) /* dot product (range [-1,1]) */
199			return( fray->rod <= -1.0 ? -1.0 :
200			fray->rod >= 1.0 ? 1.0 :
201			fray->rod );
202	greg	1.10
203	greg	1.8	if (n == 11) /* scale */
204	greg	1.12	return(funcxf.sca);
205	greg	1.8
206			if (n < 15) /* origin */
207	greg	1.12	return(funcxf.xfm[3][n-12]);
208	greg	1.8
209			if (n < 18) /* i unit vector */
210	greg	1.12	return(funcxf.xfm[0][n-15] / funcxf.sca);
211	greg	1.8
212			if (n < 21) /* j unit vector */
213	greg	1.12	return(funcxf.xfm[1][n-15] / funcxf.sca);
214	greg	1.8
215			if (n < 24) /* k unit vector */
216	greg	1.12	return(funcxf.xfm[2][n-21] / funcxf.sca);
217	greg	1.10	badchan:
218			error(USER, "illegal channel number");
219	greg	1.1	}