--- ray/src/rt/func.c	1989/04/11 13:30:29	1.4
+++ ray/src/rt/func.c	2012/06/27 16:29:26	2.27
@@ -1,104 +1,243 @@
-/* Copyright (c) 1986 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: func.c,v 2.27 2012/06/27 16:29:26 greg Exp $";
 #endif
-
 /*
  *  func.c - interface to calcomp functions.
- *
- *     4/7/86
  */
 
-#include  "ray.h"
+#include "copyright.h"
 
+#include  "ray.h"
+#include  "paths.h"
 #include  "otypes.h"
+#include  "func.h"
 
 
-typedef struct {
-	double  xfm[4][4];		/* transform matrix */
-	double  sca;			/* scalefactor */
-}  XF;
+#define  INITFILE	"rayinit.cal"
+#define  CALSUF		".cal"
+#define  LCALSUF	4
+char  REFVNAME[] = "`FILE_REFCNT";
 
-static OBJREC  *fobj;		/* current function object */
-static RAY  *fray;		/* current function ray */
-static XF  fxf;			/* current transformation */
+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, xfm, sca)		/* set channels for function call */
-OBJREC  *m;
-register RAY  *r;
-double  xfm[4][4];
-double  sca;
-{
-	extern double  l_noise3(), l_noise3a(), l_noise3b(), l_noise3c();
-	extern double  l_hermite(), l_fnoise3(), l_arg();
-	extern long  eclock;
-	static char  *initfile = "rayinit.cal";
+static char  rayinitcal[] = INITFILE;
 
-	if (initfile != NULL) {
-		loadfunc(initfile);
-		scompile(NULL, "Dx=$1;Dy=$2;Dz=$3;");
-		scompile(NULL, "Nx=$4;Ny=$5;Nz=$6;");
-		scompile(NULL, "Px=$7;Py=$8;Pz=$9;");
-		scompile(NULL, "T=$10;Rdot=$11;");
-		funset("arg", 1, l_arg);
-		funset("noise3", 3, l_noise3);
-		funset("noise3a", 3, l_noise3a);
-		funset("noise3b", 3, l_noise3b);
-		funset("noise3c", 3, l_noise3c);
-		funset("hermite", 5, l_hermite);
-		funset("fnoise3", 3, l_fnoise3);
-		initfile = NULL;
-	}
-	fobj = m;
-	fray = r;
-	fxf.sca = r->robs * sca;
-	multmat4(fxf.xfm, r->robx, xfm);
-	eclock++;		/* notify expression evaluator */
-}
+static double  l_erf(char *), l_erfc(char *), l_arg(char *);
 
 
-setfunc(m, r)				/* simplified interface to setmap */
-register OBJREC  *m;
-RAY  *r;
+void
+initfunc()	/* initialize function evaluation */
 {
-	register XF  *mxf;
+	if (!rayinitcal[0])	/* already done? */
+		return;
+	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(rayinitcal);
+	rayinitcal[0] = '\0';
+}
 
-	if ((mxf = (XF *)m->os) == NULL) {
-		register int  n = m->oargs.nsargs;
-		register char  **sa = m->oargs.sarg;
 
-		while (n > 0 && **sa != '-') {
-			n--;
-			sa++;
+MFUNC *
+getfunc(	/* get function for this modifier */
+	OBJREC  *m,
+	int  ff,
+	unsigned int  ef,
+	int  dofwd
+)
+{
+	char  sbuf[MAXSTR];
+	char  **arg;
+	MFUNC  *f;
+	int  ne, na;
+	int  i;
+					/* check to see if done already */
+	if ((f = (MFUNC *)m->os) != NULL)
+		return(f);
+	fobj = NULL; fray = NULL;
+	if (rayinitcal[0])		/* initialize on first call */
+		initfunc();
+	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");
 		}
-		mxf = (XF *)malloc(sizeof(XF));
-		if (mxf == NULL)
+	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->fxp = f->bxp = &unitxf;
+	} else {			/* get transform */
+		if ((f->bxp = (XF *)malloc(sizeof(XF))) == NULL)
 			goto memerr;
-		mxf->sca = 1.0;
-		setident4(mxf->xfm);
-		if (invxf(mxf->xfm, &mxf->sca, n, sa) != n)
+		if (invxf(f->bxp, na-i, arg+i) != na-i)
 			objerror(m, USER, "bad transform");
-		if (mxf->sca < 0.0)
-			mxf->sca = -mxf->sca;
-		m->os = (char *)mxf;
+		if (f->bxp->sca < 0.0)
+			f->bxp->sca = -f->bxp->sca;
+		if (dofwd) {			/* do both transforms */
+			if ((f->fxp = (XF *)malloc(sizeof(XF))) == NULL)
+				goto memerr;
+			xf(f->fxp, na-i, arg+i);
+			if (f->fxp->sca < 0.0)
+				f->fxp->sca = -f->fxp->sca;
+		}
 	}
-	setmap(m, r, mxf->xfm, mxf->sca);
-	return;
+	m->os = (char *)f;
+	return(f);
+toofew:
+	objerror(m, USER, "too few string arguments");
 memerr:
-	error(SYSTEM, "out of memory in setfunc");
-#undef  mxf
+	error(SYSTEM, "out of memory in getfunc");
+	return NULL; /* pro forma return */
 }
 
 
-loadfunc(fname)			/* load definition file */
-char  *fname;
+void
+freefunc(			/* free memory associated with modifier */
+	OBJREC  *m
+)
 {
-	extern char  *libpath;		/* library search path */
+	MFUNC  *f;
+	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->bxp != &unitxf)
+		free((void *)f->bxp);
+	if ((f->fxp != NULL) & (f->fxp != &unitxf))
+		free((void *)f->fxp);
+	free((void *)f);
+	m->os = NULL;
+}
+
+
+int
+setfunc(			/* set channels for function call */
+	OBJREC	*m,
+	RAY	*r
+)
+{
+	static RNUMBER	lastrno = ~0;
+	MFUNC		*f;
+					/* get function if any */
+	if ((f = (MFUNC *)m->os) == NULL)
+		objerror(m, CONSISTENCY, "setfunc called before getfunc");
+		
+	setcontext(f->ctx);		/* set evaluator context */
+					/* check to see if matrix set */
+	if ((m == fobj) & (r->rno == lastrno))
+		return(0);
+	fobj = m;
+	fray = r;
+	if (r->rox != NULL) {
+		if (f->bxp != &unitxf) {
+			funcxf.sca = r->rox->b.sca * f->bxp->sca;
+			multmat4(funcxf.xfm, r->rox->b.xfm, f->bxp->xfm);
+		} else
+			funcxf = r->rox->b;
+	} else
+		funcxf = *f->bxp;
+	lastrno = r->rno;
+	eclock++;		/* notify expression evaluator */
+	return(1);
+}
+
+
+int
+worldfunc(			/* special function context sans object */
+	char	*ctx,
+	RAY	*r
+)
+{
+	static RNUMBER	lastrno = ~0;
+
+	if (rayinitcal[0])		/* initialize on first call */
+		initfunc();
+					/* set evaluator context */
+	setcontext(ctx);
+					/* check if ray already set */
+	if ((fobj == NULL) & (r->rno == lastrno))
+		return(0);
+	fobj = NULL;
+	fray = r;
+	funcxf = unitxf;
+	lastrno = r->rno;
+	eclock++;		/* notify expression evaluator */
+	return(1);
+}
+
+
+void
+loadfunc(			/* load definition file */
+	char  *fname
+)
+{
 	char  *ffname;
 
-	if ((ffname = getpath(fname, libpath)) == NULL) {
+	if ((ffname = getpath(fname, getrlibpath(), R_OK)) == NULL) {
 		sprintf(errmsg, "cannot find function file \"%s\"", fname);
 		error(USER, errmsg);
 	}
@@ -106,12 +245,14 @@ char  *fname;
 }
 
 
-double
-l_arg()				/* return nth real argument */
+static double
+l_arg(char *nm)			/* return nth real argument */
 {
-	extern double  argument();
-	register int  n;
+	int  n;
 
+	if (fobj == NULL)
+		error(USER, "arg(n) called without a context");
+
 	n = argument(1) + .5;		/* round to integer */
 
 	if (n < 1)
@@ -125,45 +266,81 @@ l_arg()				/* return nth real argument */
 }
 
 
+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)));
+}
+
+
 double
-chanvalue(n)			/* return channel n to calcomp */
-register int  n;
+chanvalue(			/* return channel n to calcomp */
+	int  n
+)
 {
-	double  res;
-	register RAY  *r;
+	if (fray == NULL)
+		syntax("ray parameter used in constant expression");
 
-	n--;					/* for convenience */
+	if (--n < 0)
+		goto badchan;
 
-	if (n < 0 || n > 10)
-		error(USER, "illegal channel number");
+	if (n <= 2)			/* ray direction */
 
-	if (n == 9) {				/* distance */
+		return( (	fray->rdir[0]*funcxf.xfm[0][n] +
+				fray->rdir[1]*funcxf.xfm[1][n] +
+				fray->rdir[2]*funcxf.xfm[2][n]	)
+			 / funcxf.sca );
 
-		res = fray->rot;
-		for (r = fray->parent; r != NULL; r = r->parent)
-			res += r->rot;
-		res *= fxf.sca;
+	if (n <= 5)			/* surface normal */
 
-	} else if (n == 10) {			/* dot product */
+		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 );
 
-		res = fray->rod;
+	if (n <= 8)			/* intersection */
 
-	} else if (n < 3) {			/* ray direction */
-			res = (	fray->rdir[0]*fxf.xfm[0][n] +
-					fray->rdir[1]*fxf.xfm[1][n] +
-					fray->rdir[2]*fxf.xfm[2][n]	)
-				 / fxf.sca ;
-	} else if (n < 6) {			/* surface normal */
-			res = (	fray->ron[0]*fxf.xfm[0][n-3] +
-					fray->ron[1]*fxf.xfm[1][n-3] +
-					fray->ron[2]*fxf.xfm[2][n-3]	)
-				 / fxf.sca ;
-	} else {				/* intersection */
-			res =	fray->rop[0]*fxf.xfm[0][n-6] +
-					fray->rop[1]*fxf.xfm[1][n-6] +
-					fray->rop[2]*fxf.xfm[2][n-6] +
-						     fxf.xfm[3][n-6] ;
-	}
+		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] );
 
-	return(res);
+	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);
 }