src/rt/dielectric.c

#ifndef lint
static const char       RCSid[] = "$Id: dielectric.c,v 2.27 2015/05/21 13:54:59 greg Exp $";
#endif
/*
 *  dielectric.c - shading function for transparent materials.
 */

#include "copyright.h"

#include  "ray.h"
#include  "otypes.h"
#include  "rtotypes.h"
#include  "pmapmat.h"

#ifdef  DISPERSE
#include  "source.h"
static int disperse(OBJREC *m,RAY *r,FVECT vt,double tr,COLOR cet,COLOR abt);
static int lambda(OBJREC  *m, FVECT  v2, FVECT  dv, FVECT  lr);
#endif

static double mylog(double  x);


/*
 *     Explicit calculations for Fresnel's equation are performed,
 *  but only one square root computation is necessary.
 *     The index of refraction is given as a Hartmann equation
 *  with lambda0 equal to zero.  If the slope of Hartmann's
 *  equation is non-zero, the material disperses light upon
 *  refraction.  This condition is examined on rays traced to
 *  light sources.  If a ray is exiting a dielectric material, we
 *  check the sources to see if any would cause bright color to be
 *  directed to the viewer due to dispersion.  This gives colorful
 *  sparkle to crystals, etc.  (Only if DISPERSE is defined!)
 *
 *  Arguments for MAT_DIELECTRIC are:
 *      red     grn     blu     rndx    Hartmann
 *
 *  Arguments for MAT_INTERFACE are:
 *      red1    grn1    blu1    rndx1   red2    grn2    blu2    rndx2
 *
 *  The primaries are material transmission per unit length.
 *  MAT_INTERFACE uses dielectric1 for inside and dielectric2 for
 *  outside.
 */


#define  MLAMBDA        500             /* mean lambda */
#define  MAXLAMBDA      779             /* maximum lambda */
#define  MINLAMBDA      380             /* minimum lambda */

#define  MINCOS         0.997           /* minimum dot product for dispersion */

static double
mylog(          /* special log for extinction coefficients */
        double  x
)
{
        if (x < 1e-40)
                return(-100.);
        if (x >= 1.)
                return(0.);
        return(log(x));
}


int
m_dielectric(   /* color a ray which hit a dielectric interface */
        OBJREC  *m,
        RAY  *r
)
{
        double  cos1, cos2, nratio;
        COLOR  ctrans;
        COLOR  talb;
        int  hastexture;
        double  transdist=0, transtest=0;
        double  mirdist=0, mirtest=0;
        int     flatsurface;
        double  refl, trans;
        FVECT  dnorm;
        double  d1, d2;
        RAY  p;
        int  i;

        /* PMAP: skip refracted shadow or ambient ray if accounted for in
           photon map */
        if (shadowRayInPmap(r) || ambRayInPmap(r))
                return(1);
        
        if (m->oargs.nfargs != (m->otype==MAT_DIELECTRIC ? 5 : 8))
                objerror(m, USER, "bad arguments");

        raytexture(r, m->omod);                 /* get modifiers */

        if ( (hastexture = DOT(r->pert,r->pert) > FTINY*FTINY) )
                cos1 = raynormal(dnorm, r);     /* perturb normal */
        else {
                VCOPY(dnorm, r->ron);
                cos1 = r->rod;
        }
        flatsurface = r->ro != NULL && isflat(r->ro->otype) &&
                        !hastexture | (r->crtype & AMBIENT);

                                                /* index of refraction */
        if (m->otype == MAT_DIELECTRIC)
                nratio = m->oargs.farg[3] + m->oargs.farg[4]/MLAMBDA;
        else
                nratio = m->oargs.farg[3] / m->oargs.farg[7];
        
        if (cos1 < 0.0) {                       /* inside */
                hastexture = -hastexture;
                cos1 = -cos1;
                dnorm[0] = -dnorm[0];
                dnorm[1] = -dnorm[1];
                dnorm[2] = -dnorm[2];
                setcolor(r->cext, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
                                 -mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
                                 -mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
                setcolor(r->albedo, 0., 0., 0.);
                r->gecc = 0.;
                if (m->otype == MAT_INTERFACE) {
                        setcolor(ctrans,
                                -mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
                                -mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
                                -mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
                        setcolor(talb, 0., 0., 0.);
                } else {
                        copycolor(ctrans, cextinction);
                        copycolor(talb, salbedo);
                }
        } else {                                /* outside */
                nratio = 1.0 / nratio;

                setcolor(ctrans, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
                                 -mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
                                 -mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
                setcolor(talb, 0., 0., 0.);
                if (m->otype == MAT_INTERFACE) {
                        setcolor(r->cext,
                                -mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
                                -mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
                                -mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
                        setcolor(r->albedo, 0., 0., 0.);
                        r->gecc = 0.;
                }
        }

        d2 = 1.0 - nratio*nratio*(1.0 - cos1*cos1);     /* compute cos theta2 */

        if (d2 < FTINY)                 /* total reflection */

                refl = 1.0;

        else {                          /* refraction occurs */
                                        /* compute Fresnel's equations */
                cos2 = sqrt(d2);
                d1 = cos1;
                d2 = nratio*cos2;
                d1 = (d1 - d2) / (d1 + d2);
                refl = d1 * d1;

                d1 = 1.0 / cos1;
                d2 = nratio / cos2;
                d1 = (d1 - d2) / (d1 + d2);
                refl += d1 * d1;

                refl *= 0.5;
                trans = 1.0 - refl;

                trans *= nratio*nratio;         /* solid angle ratio */

                setcolor(p.rcoef, trans, trans, trans);

                if (rayorigin(&p, REFRACTED, r, p.rcoef) == 0) {

                                                /* compute refracted ray */
                        d1 = nratio*cos1 - cos2;
                        for (i = 0; i < 3; i++)
                                p.rdir[i] = nratio*r->rdir[i] + d1*dnorm[i];
                                                /* accidental reflection? */
                        if (hastexture &&
                                DOT(p.rdir,r->ron)*hastexture >= -FTINY) {
                                d1 *= (double)hastexture;
                                for (i = 0; i < 3; i++) /* ignore texture */
                                        p.rdir[i] = nratio*r->rdir[i] +
                                                        d1*r->ron[i];
                                normalize(p.rdir);      /* not exact */
                        } else
                                checknorm(p.rdir);
#ifdef  DISPERSE
                        if (m->otype != MAT_DIELECTRIC
                                        || r->rod > 0.0
                                        || r->crtype & SHADOW
                                        || !directvis
                                        || m->oargs.farg[4] == 0.0
                                        || !disperse(m, r, p.rdir,
                                                        trans, ctrans, talb))
#endif
                        {
                                copycolor(p.cext, ctrans);
                                copycolor(p.albedo, talb);
                                rayvalue(&p);
                                multcolor(p.rcol, p.rcoef);
                                addcolor(r->rcol, p.rcol);
                                                /* virtual distance */
                                if (flatsurface ||
                                        (1.-FTINY <= nratio) &
                                                (nratio <= 1.+FTINY)) {
                                        transtest = 2*bright(p.rcol);
                                        transdist = r->rot + p.rt;
                                }
                        }
                }
        }
        setcolor(p.rcoef, refl, refl, refl);

        if (!(r->crtype & SHADOW) &&
                        rayorigin(&p, REFLECTED, r, p.rcoef) == 0) {

                                        /* compute reflected ray */
                VSUM(p.rdir, r->rdir, dnorm, 2.*cos1);
                                        /* accidental penetration? */
                if (hastexture && DOT(p.rdir,r->ron)*hastexture <= FTINY)
                        VSUM(p.rdir, r->rdir, r->ron, 2.*r->rod);
                checknorm(p.rdir);
                rayvalue(&p);                   /* reflected ray value */

                multcolor(p.rcol, p.rcoef);     /* color contribution */
                addcolor(r->rcol, p.rcol);
                                                /* virtual distance */
                if (flatsurface) {
                        mirtest = 2*bright(p.rcol);
                        mirdist = r->rot + p.rt;
                }
        }
                                /* check distance to return */
        d1 = bright(r->rcol);
        if (transtest > d1)
                r->rt = transdist;
        else if (mirtest > d1)
                r->rt = mirdist;
                                /* rayvalue() computes absorption */
        return(1);
}


#ifdef  DISPERSE

static int
disperse(  /* check light sources for dispersion */
        OBJREC  *m,
        RAY  *r,
        FVECT  vt,
        double  tr,
        COLOR  cet,
        COLOR  abt
)
{
        RAY  sray;
        const RAY  *entray;
        FVECT  v1, v2, n1, n2;
        FVECT  dv, v2Xdv;
        double  v2Xdvv2Xdv;
        int  success = 0;
        SRCINDEX  si;
        FVECT  vtmp1, vtmp2;
        double  dtmp1, dtmp2;
        int  l1, l2;
        COLOR  ctmp;
        int  i;
        
        /*
         *     This routine computes dispersion to the first order using
         *  the following assumptions:
         *
         *      1) The dependency of the index of refraction on wavelength
         *              is approximated by Hartmann's equation with lambda0
         *              equal to zero.
         *      2) The entry and exit locations are constant with respect
         *              to dispersion.
         *
         *     The second assumption permits us to model dispersion without
         *  having to sample refracted directions.  We assume that the
         *  geometry inside the material is constant, and concern ourselves
         *  only with the relationship between the entering and exiting ray.
         *  We compute the first derivatives of the entering and exiting
         *  refraction with respect to the index of refraction.  This 
         *  is then used in a first order Taylor series to determine the
         *  index of refraction necessary to send the exiting ray to each
         *  light source.
         *     If an exiting ray hits a light source within the refraction
         *  boundaries, we sum all the frequencies over the disc of the
         *  light source to determine the resulting color.  A smaller light
         *  source will therefore exhibit a sharper spectrum.
         */

        if (!(r->crtype & REFRACTED)) {         /* ray started in material */
                VCOPY(v1, r->rdir);
                n1[0] = -r->rdir[0]; n1[1] = -r->rdir[1]; n1[2] = -r->rdir[2];
        } else {
                                                /* find entry point */
                for (entray = r; entray->rtype != REFRACTED;
                                entray = entray->parent)
                        ;
                entray = entray->parent;
                if (entray->crtype & REFRACTED)         /* too difficult */
                        return(0);
                VCOPY(v1, entray->rdir);
                VCOPY(n1, entray->ron);
        }
        VCOPY(v2, vt);                  /* exiting ray */
        VCOPY(n2, r->ron);

                                        /* first order dispersion approx. */
        dtmp1 = 1./DOT(n1, v1);
        dtmp2 = 1./DOT(n2, v2);
        for (i = 0; i < 3; i++)
                dv[i] = v1[i] + v2[i] - n1[i]*dtmp1 - n2[i]*dtmp2;
                
        if (DOT(dv, dv) <= FTINY)       /* null effect */
                return(0);
                                        /* compute plane normal */
        fcross(v2Xdv, v2, dv);
        v2Xdvv2Xdv = DOT(v2Xdv, v2Xdv);

                                        /* check sources */
        initsrcindex(&si);
        while (srcray(&sray, r, &si)) {

                if (DOT(sray.rdir, v2) < MINCOS)
                        continue;                       /* bad source */
                                                /* adjust source ray */

                dtmp1 = DOT(v2Xdv, sray.rdir) / v2Xdvv2Xdv;
                sray.rdir[0] -= dtmp1 * v2Xdv[0];
                sray.rdir[1] -= dtmp1 * v2Xdv[1];
                sray.rdir[2] -= dtmp1 * v2Xdv[2];

                l1 = lambda(m, v2, dv, sray.rdir);      /* mean lambda */

                if (l1 > MAXLAMBDA || l1 < MINLAMBDA)   /* not visible */
                        continue;
                                                /* trace source ray */
                copycolor(sray.cext, cet);
                copycolor(sray.albedo, abt);
                normalize(sray.rdir);
                rayvalue(&sray);
                if (bright(sray.rcol) <= FTINY) /* missed it */
                        continue;
                
                /*
                 *      Compute spectral sum over diameter of source.
                 *  First find directions for rays going to opposite
                 *  sides of source, then compute wavelengths for each.
                 */
                 
                fcross(vtmp1, v2Xdv, sray.rdir);
                dtmp1 = sqrt(si.dom  / v2Xdvv2Xdv / PI);

                                                        /* compute first ray */
                VSUM(vtmp2, sray.rdir, vtmp1, dtmp1);

                l1 = lambda(m, v2, dv, vtmp2);          /* first lambda */
                if (l1 < 0)
                        continue;
                                                        /* compute second ray */
                VSUM(vtmp2, sray.rdir, vtmp1, -dtmp1);

                l2 = lambda(m, v2, dv, vtmp2);          /* second lambda */
                if (l2 < 0)
                        continue;
                                        /* compute color from spectrum */
                if (l1 < l2)
                        spec_rgb(ctmp, l1, l2);
                else
                        spec_rgb(ctmp, l2, l1);
                multcolor(ctmp, sray.rcol);
                scalecolor(ctmp, tr);
                addcolor(r->rcol, ctmp);
                success++;
        }
        return(success);
}


static int
lambda(                 /* compute lambda for material */
        OBJREC  *m,
        FVECT  v2,
        FVECT  dv,
        FVECT  lr
)
{
        FVECT  lrXdv, v2Xlr;
        double  dtmp, denom;
        int  i;

        fcross(lrXdv, lr, dv);
        for (i = 0; i < 3; i++) 
                if ((lrXdv[i] > FTINY) | (lrXdv[i] < -FTINY))
                        break;
        if (i >= 3)
                return(-1);

        fcross(v2Xlr, v2, lr);

        dtmp = m->oargs.farg[4] / MLAMBDA;
        denom = dtmp + v2Xlr[i]/lrXdv[i] * (m->oargs.farg[3] + dtmp);

        if (denom < FTINY)
                return(-1);

        return(m->oargs.farg[4] / denom);
}

#endif  /* DISPERSE */
Revision:	2.28
Committed:	Wed Oct 28 15:45:58 2015 UTC (8 years, 6 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R2, rad5R1
Changes since 2.27:	+2 -2 lines
Log Message:	Added back ambient ray testing for photon map, which is needed by rcontrib
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.28	static const char RCSid[] = "$Id: dielectric.c,v 2.27 2015/05/21 13:54:59 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* dielectric.c - shading function for transparent materials.
6	greg	2.15	*/
7
8	greg	2.16	#include "copyright.h"
9	greg	1.1
10			#include "ray.h"
11			#include "otypes.h"
12	schorsch	2.18	#include "rtotypes.h"
13	greg	2.25	#include "pmapmat.h"
14	greg	1.1
15			#ifdef DISPERSE
16			#include "source.h"
17	schorsch	2.18	static int disperse(OBJREC m,RAY r,FVECT vt,double tr,COLOR cet,COLOR abt);
18			static int lambda(OBJREC *m, FVECT v2, FVECT dv, FVECT lr);
19	greg	1.1	#endif
20
21	schorsch	2.18	static double mylog(double x);
22
23
24	greg	1.1	/*
25			* Explicit calculations for Fresnel's equation are performed,
26			* but only one square root computation is necessary.
27			* The index of refraction is given as a Hartmann equation
28			* with lambda0 equal to zero. If the slope of Hartmann's
29			* equation is non-zero, the material disperses light upon
30			* refraction. This condition is examined on rays traced to
31			* light sources. If a ray is exiting a dielectric material, we
32			* check the sources to see if any would cause bright color to be
33			* directed to the viewer due to dispersion. This gives colorful
34			* sparkle to crystals, etc. (Only if DISPERSE is defined!)
35			*
36			* Arguments for MAT_DIELECTRIC are:
37			* red grn blu rndx Hartmann
38			*
39			* Arguments for MAT_INTERFACE are:
40			* red1 grn1 blu1 rndx1 red2 grn2 blu2 rndx2
41			*
42			* The primaries are material transmission per unit length.
43			* MAT_INTERFACE uses dielectric1 for inside and dielectric2 for
44			* outside.
45			*/
46
47
48			#define MLAMBDA 500 /* mean lambda */
49			#define MAXLAMBDA 779 /* maximum lambda */
50			#define MINLAMBDA 380 /* minimum lambda */
51
52			#define MINCOS 0.997 /* minimum dot product for dispersion */
53
54	greg	2.27	static double
55	schorsch	2.18	mylog( /* special log for extinction coefficients */
56			double x
57			)
58	greg	2.10	{
59			if (x < 1e-40)
60			return(-100.);
61			if (x >= 1.)
62			return(0.);
63			return(log(x));
64			}
65
66
67	greg	2.24	int
68	schorsch	2.18	m_dielectric( /* color a ray which hit a dielectric interface */
69			OBJREC *m,
70	greg	2.24	RAY *r
71	schorsch	2.18	)
72	greg	1.1	{
73			double cos1, cos2, nratio;
74	greg	2.9	COLOR ctrans;
75	greg	2.11	COLOR talb;
76	gwlarson	2.14	int hastexture;
77	greg	2.23	double transdist=0, transtest=0;
78			double mirdist=0, mirtest=0;
79	greg	2.19	int flatsurface;
80	greg	1.5	double refl, trans;
81	greg	1.1	FVECT dnorm;
82			double d1, d2;
83			RAY p;
84	greg	2.24	int i;
85	greg	1.1
86	rschregle	2.26	/* PMAP: skip refracted shadow or ambient ray if accounted for in
87			photon map */
88	greg	2.28	if (shadowRayInPmap(r) \|\| ambRayInPmap(r))
89	greg	2.25	return(1);
90
91	greg	1.1	if (m->oargs.nfargs != (m->otype==MAT_DIELECTRIC ? 5 : 8))
92			objerror(m, USER, "bad arguments");
93
94			raytexture(r, m->omod); /* get modifiers */
95
96	schorsch	2.17	if ( (hastexture = DOT(r->pert,r->pert) > FTINY*FTINY) )
97	gwlarson	2.14	cos1 = raynormal(dnorm, r); /* perturb normal */
98			else {
99			VCOPY(dnorm, r->ron);
100			cos1 = r->rod;
101			}
102	greg	2.24	flatsurface = r->ro != NULL && isflat(r->ro->otype) &&
103			!hastexture \| (r->crtype & AMBIENT);
104	greg	2.19
105	greg	1.1	/* index of refraction */
106			if (m->otype == MAT_DIELECTRIC)
107			nratio = m->oargs.farg[3] + m->oargs.farg[4]/MLAMBDA;
108			else
109			nratio = m->oargs.farg[3] / m->oargs.farg[7];
110
111			if (cos1 < 0.0) { /* inside */
112	gwlarson	2.14	hastexture = -hastexture;
113	greg	1.1	cos1 = -cos1;
114			dnorm[0] = -dnorm[0];
115			dnorm[1] = -dnorm[1];
116			dnorm[2] = -dnorm[2];
117	greg	2.10	setcolor(r->cext, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
118			-mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
119			-mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
120	greg	2.11	setcolor(r->albedo, 0., 0., 0.);
121	greg	2.9	r->gecc = 0.;
122			if (m->otype == MAT_INTERFACE) {
123			setcolor(ctrans,
124	greg	2.10	-mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
125			-mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
126			-mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
127	greg	2.11	setcolor(talb, 0., 0., 0.);
128	greg	2.9	} else {
129			copycolor(ctrans, cextinction);
130	greg	2.11	copycolor(talb, salbedo);
131	greg	2.9	}
132	greg	1.1	} else { /* outside */
133			nratio = 1.0 / nratio;
134	greg	2.9
135	greg	2.10	setcolor(ctrans, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
136			-mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
137			-mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
138	greg	2.11	setcolor(talb, 0., 0., 0.);
139	greg	2.9	if (m->otype == MAT_INTERFACE) {
140			setcolor(r->cext,
141	greg	2.10	-mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
142			-mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
143			-mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
144	greg	2.11	setcolor(r->albedo, 0., 0., 0.);
145	greg	2.9	r->gecc = 0.;
146			}
147	greg	1.1	}
148
149			d2 = 1.0 - nrationratio(1.0 - cos1cos1); / compute cos theta2 */
150
151			if (d2 < FTINY) /* total reflection */
152
153			refl = 1.0;
154
155			else { /* refraction occurs */
156			/* compute Fresnel's equations */
157			cos2 = sqrt(d2);
158			d1 = cos1;
159			d2 = nratio*cos2;
160			d1 = (d1 - d2) / (d1 + d2);
161			refl = d1 * d1;
162
163			d1 = 1.0 / cos1;
164			d2 = nratio / cos2;
165			d1 = (d1 - d2) / (d1 + d2);
166			refl += d1 * d1;
167
168	greg	2.9	refl *= 0.5;
169	greg	1.1	trans = 1.0 - refl;
170	greg	2.15
171			trans = nrationratio; /* solid angle ratio */
172	greg	1.1
173	greg	2.20	setcolor(p.rcoef, trans, trans, trans);
174
175			if (rayorigin(&p, REFRACTED, r, p.rcoef) == 0) {
176	greg	1.1
177			/* compute refracted ray */
178			d1 = nratio*cos1 - cos2;
179			for (i = 0; i < 3; i++)
180			p.rdir[i] = nratior->rdir[i] + d1dnorm[i];
181	gwlarson	2.14	/* accidental reflection? */
182			if (hastexture &&
183			DOT(p.rdir,r->ron)*hastexture >= -FTINY) {
184			d1 *= (double)hastexture;
185			for (i = 0; i < 3; i++) /* ignore texture */
186			p.rdir[i] = nratio*r->rdir[i] +
187			d1*r->ron[i];
188			normalize(p.rdir); /* not exact */
189	greg	2.21	} else
190			checknorm(p.rdir);
191	greg	1.1	#ifdef DISPERSE
192			if (m->otype != MAT_DIELECTRIC
193			\|\| r->rod > 0.0
194			\|\| r->crtype & SHADOW
195	greg	2.3	\|\| !directvis
196	greg	1.1	\|\| m->oargs.farg[4] == 0.0
197	greg	2.12	\|\| !disperse(m, r, p.rdir,
198			trans, ctrans, talb))
199	greg	1.1	#endif
200			{
201	greg	2.9	copycolor(p.cext, ctrans);
202	greg	2.11	copycolor(p.albedo, talb);
203	greg	1.1	rayvalue(&p);
204	greg	2.20	multcolor(p.rcol, p.rcoef);
205	greg	1.1	addcolor(r->rcol, p.rcol);
206	greg	2.19	/* virtual distance */
207			if (flatsurface \|\|
208	greg	2.24	(1.-FTINY <= nratio) &
209			(nratio <= 1.+FTINY)) {
210	greg	2.19	transtest = 2*bright(p.rcol);
211			transdist = r->rot + p.rt;
212			}
213	greg	1.1	}
214			}
215			}
216	greg	2.20	setcolor(p.rcoef, refl, refl, refl);
217
218	greg	1.1	if (!(r->crtype & SHADOW) &&
219	greg	2.20	rayorigin(&p, REFLECTED, r, p.rcoef) == 0) {
220	greg	1.1
221			/* compute reflected ray */
222	greg	2.22	VSUM(p.rdir, r->rdir, dnorm, 2.*cos1);
223	gwlarson	2.14	/* accidental penetration? */
224			if (hastexture && DOT(p.rdir,r->ron)*hastexture <= FTINY)
225	greg	2.22	VSUM(p.rdir, r->rdir, r->ron, 2.*r->rod);
226	greg	2.21	checknorm(p.rdir);
227	greg	1.1	rayvalue(&p); /* reflected ray value */
228
229	greg	2.20	multcolor(p.rcol, p.rcoef); /* color contribution */
230	greg	1.1	addcolor(r->rcol, p.rcol);
231	greg	2.19	/* virtual distance */
232			if (flatsurface) {
233			mirtest = 2*bright(p.rcol);
234			mirdist = r->rot + p.rt;
235			}
236	greg	1.1	}
237	greg	2.19	/* check distance to return */
238			d1 = bright(r->rcol);
239			if (transtest > d1)
240			r->rt = transdist;
241			else if (mirtest > d1)
242			r->rt = mirdist;
243	greg	2.9	/* rayvalue() computes absorption */
244	greg	2.7	return(1);
245	greg	1.1	}
246
247
248			#ifdef DISPERSE
249
250	schorsch	2.18	static int
251			disperse( /* check light sources for dispersion */
252			OBJREC *m,
253			RAY *r,
254			FVECT vt,
255			double tr,
256			COLOR cet,
257			COLOR abt
258			)
259	greg	1.1	{
260	greg	2.20	RAY sray;
261			const RAY *entray;
262	greg	1.1	FVECT v1, v2, n1, n2;
263			FVECT dv, v2Xdv;
264			double v2Xdvv2Xdv;
265	greg	1.7	int success = 0;
266			SRCINDEX si;
267	greg	1.1	FVECT vtmp1, vtmp2;
268			double dtmp1, dtmp2;
269			int l1, l2;
270			COLOR ctmp;
271			int i;
272
273			/*
274			* This routine computes dispersion to the first order using
275			* the following assumptions:
276			*
277			* 1) The dependency of the index of refraction on wavelength
278			* is approximated by Hartmann's equation with lambda0
279			* equal to zero.
280			* 2) The entry and exit locations are constant with respect
281			* to dispersion.
282			*
283			* The second assumption permits us to model dispersion without
284			* having to sample refracted directions. We assume that the
285			* geometry inside the material is constant, and concern ourselves
286			* only with the relationship between the entering and exiting ray.
287			* We compute the first derivatives of the entering and exiting
288			* refraction with respect to the index of refraction. This
289			* is then used in a first order Taylor series to determine the
290			* index of refraction necessary to send the exiting ray to each
291			* light source.
292			* If an exiting ray hits a light source within the refraction
293			* boundaries, we sum all the frequencies over the disc of the
294			* light source to determine the resulting color. A smaller light
295			* source will therefore exhibit a sharper spectrum.
296			*/
297
298			if (!(r->crtype & REFRACTED)) { /* ray started in material */
299			VCOPY(v1, r->rdir);
300			n1[0] = -r->rdir[0]; n1[1] = -r->rdir[1]; n1[2] = -r->rdir[2];
301			} else {
302			/* find entry point */
303			for (entray = r; entray->rtype != REFRACTED;
304			entray = entray->parent)
305			;
306			entray = entray->parent;
307			if (entray->crtype & REFRACTED) /* too difficult */
308			return(0);
309			VCOPY(v1, entray->rdir);
310			VCOPY(n1, entray->ron);
311			}
312			VCOPY(v2, vt); /* exiting ray */
313			VCOPY(n2, r->ron);
314
315			/* first order dispersion approx. */
316	greg	2.22	dtmp1 = 1./DOT(n1, v1);
317			dtmp2 = 1./DOT(n2, v2);
318	greg	1.1	for (i = 0; i < 3; i++)
319	greg	2.22	dv[i] = v1[i] + v2[i] - n1[i]dtmp1 - n2[i]dtmp2;
320	greg	1.1
321			if (DOT(dv, dv) <= FTINY) /* null effect */
322			return(0);
323			/* compute plane normal */
324			fcross(v2Xdv, v2, dv);
325			v2Xdvv2Xdv = DOT(v2Xdv, v2Xdv);
326
327			/* check sources */
328	greg	1.7	initsrcindex(&si);
329			while (srcray(&sray, r, &si)) {
330	greg	1.1
331	greg	1.7	if (DOT(sray.rdir, v2) < MINCOS)
332	greg	1.1	continue; /* bad source */
333			/* adjust source ray */
334
335			dtmp1 = DOT(v2Xdv, sray.rdir) / v2Xdvv2Xdv;
336			sray.rdir[0] -= dtmp1 * v2Xdv[0];
337			sray.rdir[1] -= dtmp1 * v2Xdv[1];
338			sray.rdir[2] -= dtmp1 * v2Xdv[2];
339
340			l1 = lambda(m, v2, dv, sray.rdir); /* mean lambda */
341
342			if (l1 > MAXLAMBDA \|\| l1 < MINLAMBDA) /* not visible */
343			continue;
344			/* trace source ray */
345	greg	2.12	copycolor(sray.cext, cet);
346			copycolor(sray.albedo, abt);
347	greg	1.1	normalize(sray.rdir);
348			rayvalue(&sray);
349	greg	1.2	if (bright(sray.rcol) <= FTINY) /* missed it */
350	greg	1.1	continue;
351
352			/*
353			* Compute spectral sum over diameter of source.
354			* First find directions for rays going to opposite
355			* sides of source, then compute wavelengths for each.
356			*/
357
358			fcross(vtmp1, v2Xdv, sray.rdir);
359	greg	1.7	dtmp1 = sqrt(si.dom / v2Xdvv2Xdv / PI);
360	greg	1.1
361			/* compute first ray */
362	greg	2.22	VSUM(vtmp2, sray.rdir, vtmp1, dtmp1);
363	greg	1.1
364			l1 = lambda(m, v2, dv, vtmp2); /* first lambda */
365			if (l1 < 0)
366			continue;
367			/* compute second ray */
368	greg	2.22	VSUM(vtmp2, sray.rdir, vtmp1, -dtmp1);
369	greg	1.1
370			l2 = lambda(m, v2, dv, vtmp2); /* second lambda */
371			if (l2 < 0)
372			continue;
373			/* compute color from spectrum */
374			if (l1 < l2)
375			spec_rgb(ctmp, l1, l2);
376			else
377			spec_rgb(ctmp, l2, l1);
378			multcolor(ctmp, sray.rcol);
379			scalecolor(ctmp, tr);
380			addcolor(r->rcol, ctmp);
381			success++;
382			}
383			return(success);
384			}
385
386
387			static int
388	schorsch	2.18	lambda( /* compute lambda for material */
389	greg	2.24	OBJREC *m,
390	schorsch	2.18	FVECT v2,
391			FVECT dv,
392			FVECT lr
393			)
394	greg	1.1	{
395			FVECT lrXdv, v2Xlr;
396			double dtmp, denom;
397			int i;
398
399			fcross(lrXdv, lr, dv);
400			for (i = 0; i < 3; i++)
401	greg	2.22	if ((lrXdv[i] > FTINY) \| (lrXdv[i] < -FTINY))
402	greg	1.1	break;
403			if (i >= 3)
404			return(-1);
405
406			fcross(v2Xlr, v2, lr);
407
408			dtmp = m->oargs.farg[4] / MLAMBDA;
409			denom = dtmp + v2Xlr[i]/lrXdv[i] * (m->oargs.farg[3] + dtmp);
410
411			if (denom < FTINY)
412			return(-1);
413
414			return(m->oargs.farg[4] / denom);
415			}
416
417			#endif /* DISPERSE */