src/rt/dielectric.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#endif
/*
 *  dielectric.c - shading function for transparent materials.
 */

/* ====================================================================
 * The Radiance Software License, Version 1.0
 *
 * Copyright (c) 1990 - 2002 The Regents of the University of California,
 * through Lawrence Berkeley National Laboratory.   All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *         notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in
 *       the documentation and/or other materials provided with the
 *       distribution.
 *
 * 3. The end-user documentation included with the redistribution,
 *           if any, must include the following acknowledgment:
 *             "This product includes Radiance software
 *                 (http://radsite.lbl.gov/)
 *                 developed by the Lawrence Berkeley National Laboratory
 *               (http://www.lbl.gov/)."
 *       Alternately, this acknowledgment may appear in the software itself,
 *       if and wherever such third-party acknowledgments normally appear.
 *
 * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
 *       and "The Regents of the University of California" must
 *       not be used to endorse or promote products derived from this
 *       software without prior written permission. For written
 *       permission, please contact [email protected].
 *
 * 5. Products derived from this software may not be called "Radiance",
 *       nor may "Radiance" appear in their name, without prior written
 *       permission of Lawrence Berkeley National Laboratory.
 *
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED.   IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
 * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
 * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 * ====================================================================
 *
 * This software consists of voluntary contributions made by many
 * individuals on behalf of Lawrence Berkeley National Laboratory.   For more
 * information on Lawrence Berkeley National Laboratory, please see
 * <http://www.lbl.gov/>.
 */

#include  "ray.h"

#include  "otypes.h"

#ifdef  DISPERSE
#include  "source.h"
static  disperse();
static int  lambda();
#endif

/*
 *     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(x)                /* special log for extinction coefficients */
double  x;
{
        if (x < 1e-40)
                return(-100.);
        if (x >= 1.)
                return(0.);
        return(log(x));
}


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

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

                if (rayorigin(&p, r, REFRACTED, trans) == 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 */
                        }
#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);
                                scalecolor(p.rcol, trans);
                                addcolor(r->rcol, p.rcol);
                                if (nratio >= 1.0-FTINY && nratio <= 1.0+FTINY)
                                        r->rt = r->rot + p.rt;
                        }
                }
        }
                
        if (!(r->crtype & SHADOW) &&
                        rayorigin(&p, r, REFLECTED, refl) == 0) {

                                        /* compute reflected ray */
                for (i = 0; i < 3; i++)
                        p.rdir[i] = r->rdir[i] + 2.0*cos1*dnorm[i];
                                        /* accidental penetration? */
                if (hastexture && DOT(p.rdir,r->ron)*hastexture <= FTINY)
                        for (i = 0; i < 3; i++)         /* ignore texture */
                                p.rdir[i] = r->rdir[i] + 2.0*r->rod*r->ron[i];

                rayvalue(&p);                   /* reflected ray value */

                scalecolor(p.rcol, refl);       /* color contribution */
                addcolor(r->rcol, p.rcol);
        }
                                /* rayvalue() computes absorption */
        return(1);
}


#ifdef  DISPERSE

static
disperse(m, r, vt, tr, cet, abt)  /* check light sources for dispersion */
OBJREC  *m;
RAY  *r;
FVECT  vt;
double  tr;
COLOR  cet, abt;
{
        RAY  sray, *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 = DOT(n1, v1);
        dtmp2 = 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 */
                for (i = 0; i < 3; i++)
                        vtmp2[i] = sray.rdir[i] + dtmp1*vtmp1[i];

                l1 = lambda(m, v2, dv, vtmp2);          /* first lambda */
                if (l1 < 0)
                        continue;
                                                        /* compute second ray */
                for (i = 0; i < 3; i++)
                        vtmp2[i] = sray.rdir[i] - dtmp1*vtmp1[i];

                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(m, v2, dv, lr)                   /* compute lambda for material */
register OBJREC  *m;
FVECT  v2, dv, 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.15
Committed:	Sat Feb 22 02:07:28 2003 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.14:	+59 -8 lines
Log Message:	Changes and check-in for 3.5 release Includes new source files and modifications not recorded for many years See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.15	static const char RCSid[] = "$Id$";
3	greg	1.1	#endif
4			/*
5			* dielectric.c - shading function for transparent materials.
6	greg	2.15	*/
7
8			/* ====================================================================
9			* The Radiance Software License, Version 1.0
10			*
11			* Copyright (c) 1990 - 2002 The Regents of the University of California,
12			* through Lawrence Berkeley National Laboratory. All rights reserved.
13			*
14			* Redistribution and use in source and binary forms, with or without
15			* modification, are permitted provided that the following conditions
16			* are met:
17			*
18			* 1. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20	greg	1.1	*
21	greg	2.15	* 2. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in
23			* the documentation and/or other materials provided with the
24			* distribution.
25			*
26			* 3. The end-user documentation included with the redistribution,
27			* if any, must include the following acknowledgment:
28			* "This product includes Radiance software
29			* (http://radsite.lbl.gov/)
30			* developed by the Lawrence Berkeley National Laboratory
31			* (http://www.lbl.gov/)."
32			* Alternately, this acknowledgment may appear in the software itself,
33			* if and wherever such third-party acknowledgments normally appear.
34			*
35			* 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
36			* and "The Regents of the University of California" must
37			* not be used to endorse or promote products derived from this
38			* software without prior written permission. For written
39			* permission, please contact [email protected].
40			*
41			* 5. Products derived from this software may not be called "Radiance",
42			* nor may "Radiance" appear in their name, without prior written
43			* permission of Lawrence Berkeley National Laboratory.
44			*
45			* THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
46			* WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
47			* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
48			* DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
49			* ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
50			* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
51			* LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
52			* USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
53			* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
54			* OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
55			* OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
56			* SUCH DAMAGE.
57			* ====================================================================
58			*
59			* This software consists of voluntary contributions made by many
60			* individuals on behalf of Lawrence Berkeley National Laboratory. For more
61			* information on Lawrence Berkeley National Laboratory, please see
62			* <http://www.lbl.gov/>.
63	greg	1.1	*/
64
65			#include "ray.h"
66
67			#include "otypes.h"
68
69			#ifdef DISPERSE
70			#include "source.h"
71	greg	2.5	static disperse();
72	greg	2.6	static int lambda();
73	greg	1.1	#endif
74
75			/*
76			* Explicit calculations for Fresnel's equation are performed,
77			* but only one square root computation is necessary.
78			* The index of refraction is given as a Hartmann equation
79			* with lambda0 equal to zero. If the slope of Hartmann's
80			* equation is non-zero, the material disperses light upon
81			* refraction. This condition is examined on rays traced to
82			* light sources. If a ray is exiting a dielectric material, we
83			* check the sources to see if any would cause bright color to be
84			* directed to the viewer due to dispersion. This gives colorful
85			* sparkle to crystals, etc. (Only if DISPERSE is defined!)
86			*
87			* Arguments for MAT_DIELECTRIC are:
88			* red grn blu rndx Hartmann
89			*
90			* Arguments for MAT_INTERFACE are:
91			* red1 grn1 blu1 rndx1 red2 grn2 blu2 rndx2
92			*
93			* The primaries are material transmission per unit length.
94			* MAT_INTERFACE uses dielectric1 for inside and dielectric2 for
95			* outside.
96			*/
97
98
99			#define MLAMBDA 500 /* mean lambda */
100			#define MAXLAMBDA 779 /* maximum lambda */
101			#define MINLAMBDA 380 /* minimum lambda */
102
103			#define MINCOS 0.997 /* minimum dot product for dispersion */
104
105	greg	2.9
106	greg	2.10	static double
107			mylog(x) /* special log for extinction coefficients */
108			double x;
109			{
110			if (x < 1e-40)
111			return(-100.);
112			if (x >= 1.)
113			return(0.);
114			return(log(x));
115			}
116
117
118	greg	2.9	m_dielectric(m, r) /* color a ray which hit a dielectric interface */
119	greg	1.1	OBJREC *m;
120			register RAY *r;
121			{
122			double cos1, cos2, nratio;
123	greg	2.9	COLOR ctrans;
124	greg	2.11	COLOR talb;
125	gwlarson	2.14	int hastexture;
126	greg	1.5	double refl, trans;
127	greg	1.1	FVECT dnorm;
128			double d1, d2;
129			RAY p;
130			register int i;
131
132			if (m->oargs.nfargs != (m->otype==MAT_DIELECTRIC ? 5 : 8))
133			objerror(m, USER, "bad arguments");
134
135			raytexture(r, m->omod); /* get modifiers */
136
137	gwlarson	2.14	if (hastexture = DOT(r->pert,r->pert) > FTINY*FTINY)
138			cos1 = raynormal(dnorm, r); /* perturb normal */
139			else {
140			VCOPY(dnorm, r->ron);
141			cos1 = r->rod;
142			}
143	greg	1.1	/* index of refraction */
144			if (m->otype == MAT_DIELECTRIC)
145			nratio = m->oargs.farg[3] + m->oargs.farg[4]/MLAMBDA;
146			else
147			nratio = m->oargs.farg[3] / m->oargs.farg[7];
148
149			if (cos1 < 0.0) { /* inside */
150	gwlarson	2.14	hastexture = -hastexture;
151	greg	1.1	cos1 = -cos1;
152			dnorm[0] = -dnorm[0];
153			dnorm[1] = -dnorm[1];
154			dnorm[2] = -dnorm[2];
155	greg	2.10	setcolor(r->cext, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
156			-mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
157			-mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
158	greg	2.11	setcolor(r->albedo, 0., 0., 0.);
159	greg	2.9	r->gecc = 0.;
160			if (m->otype == MAT_INTERFACE) {
161			setcolor(ctrans,
162	greg	2.10	-mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
163			-mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
164			-mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
165	greg	2.11	setcolor(talb, 0., 0., 0.);
166	greg	2.9	} else {
167			copycolor(ctrans, cextinction);
168	greg	2.11	copycolor(talb, salbedo);
169	greg	2.9	}
170	greg	1.1	} else { /* outside */
171			nratio = 1.0 / nratio;
172	greg	2.9
173	greg	2.10	setcolor(ctrans, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
174			-mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
175			-mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
176	greg	2.11	setcolor(talb, 0., 0., 0.);
177	greg	2.9	if (m->otype == MAT_INTERFACE) {
178			setcolor(r->cext,
179	greg	2.10	-mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
180			-mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
181			-mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
182	greg	2.11	setcolor(r->albedo, 0., 0., 0.);
183	greg	2.9	r->gecc = 0.;
184			}
185	greg	1.1	}
186
187			d2 = 1.0 - nrationratio(1.0 - cos1cos1); / compute cos theta2 */
188
189			if (d2 < FTINY) /* total reflection */
190
191			refl = 1.0;
192
193			else { /* refraction occurs */
194			/* compute Fresnel's equations */
195			cos2 = sqrt(d2);
196			d1 = cos1;
197			d2 = nratio*cos2;
198			d1 = (d1 - d2) / (d1 + d2);
199			refl = d1 * d1;
200
201			d1 = 1.0 / cos1;
202			d2 = nratio / cos2;
203			d1 = (d1 - d2) / (d1 + d2);
204			refl += d1 * d1;
205
206	greg	2.9	refl *= 0.5;
207	greg	1.1	trans = 1.0 - refl;
208	greg	2.15
209			trans = nrationratio; /* solid angle ratio */
210	greg	1.1
211	gwlarson	2.13	if (rayorigin(&p, r, REFRACTED, trans) == 0) {
212	greg	1.1
213			/* compute refracted ray */
214			d1 = nratio*cos1 - cos2;
215			for (i = 0; i < 3; i++)
216			p.rdir[i] = nratior->rdir[i] + d1dnorm[i];
217	gwlarson	2.14	/* accidental reflection? */
218			if (hastexture &&
219			DOT(p.rdir,r->ron)*hastexture >= -FTINY) {
220			d1 *= (double)hastexture;
221			for (i = 0; i < 3; i++) /* ignore texture */
222			p.rdir[i] = nratio*r->rdir[i] +
223			d1*r->ron[i];
224			normalize(p.rdir); /* not exact */
225			}
226	greg	1.1	#ifdef DISPERSE
227			if (m->otype != MAT_DIELECTRIC
228			\|\| r->rod > 0.0
229			\|\| r->crtype & SHADOW
230	greg	2.3	\|\| !directvis
231	greg	1.1	\|\| m->oargs.farg[4] == 0.0
232	greg	2.12	\|\| !disperse(m, r, p.rdir,
233			trans, ctrans, talb))
234	greg	1.1	#endif
235			{
236	greg	2.9	copycolor(p.cext, ctrans);
237	greg	2.11	copycolor(p.albedo, talb);
238	greg	1.1	rayvalue(&p);
239			scalecolor(p.rcol, trans);
240			addcolor(r->rcol, p.rcol);
241	greg	2.4	if (nratio >= 1.0-FTINY && nratio <= 1.0+FTINY)
242			r->rt = r->rot + p.rt;
243	greg	1.1	}
244			}
245			}
246
247			if (!(r->crtype & SHADOW) &&
248	gwlarson	2.13	rayorigin(&p, r, REFLECTED, refl) == 0) {
249	greg	1.1
250			/* compute reflected ray */
251			for (i = 0; i < 3; i++)
252			p.rdir[i] = r->rdir[i] + 2.0cos1dnorm[i];
253	gwlarson	2.14	/* accidental penetration? */
254			if (hastexture && DOT(p.rdir,r->ron)*hastexture <= FTINY)
255			for (i = 0; i < 3; i++) /* ignore texture */
256			p.rdir[i] = r->rdir[i] + 2.0r->rodr->ron[i];
257	greg	1.1
258			rayvalue(&p); /* reflected ray value */
259
260			scalecolor(p.rcol, refl); /* color contribution */
261			addcolor(r->rcol, p.rcol);
262			}
263	greg	2.9	/* rayvalue() computes absorption */
264	greg	2.7	return(1);
265	greg	1.1	}
266
267
268			#ifdef DISPERSE
269
270			static
271	greg	2.12	disperse(m, r, vt, tr, cet, abt) /* check light sources for dispersion */
272	greg	1.1	OBJREC *m;
273			RAY *r;
274			FVECT vt;
275			double tr;
276	greg	2.12	COLOR cet, abt;
277	greg	1.1	{
278			RAY sray, *entray;
279			FVECT v1, v2, n1, n2;
280			FVECT dv, v2Xdv;
281			double v2Xdvv2Xdv;
282	greg	1.7	int success = 0;
283			SRCINDEX si;
284	greg	1.1	FVECT vtmp1, vtmp2;
285			double dtmp1, dtmp2;
286			int l1, l2;
287			COLOR ctmp;
288			int i;
289
290			/*
291			* This routine computes dispersion to the first order using
292			* the following assumptions:
293			*
294			* 1) The dependency of the index of refraction on wavelength
295			* is approximated by Hartmann's equation with lambda0
296			* equal to zero.
297			* 2) The entry and exit locations are constant with respect
298			* to dispersion.
299			*
300			* The second assumption permits us to model dispersion without
301			* having to sample refracted directions. We assume that the
302			* geometry inside the material is constant, and concern ourselves
303			* only with the relationship between the entering and exiting ray.
304			* We compute the first derivatives of the entering and exiting
305			* refraction with respect to the index of refraction. This
306			* is then used in a first order Taylor series to determine the
307			* index of refraction necessary to send the exiting ray to each
308			* light source.
309			* If an exiting ray hits a light source within the refraction
310			* boundaries, we sum all the frequencies over the disc of the
311			* light source to determine the resulting color. A smaller light
312			* source will therefore exhibit a sharper spectrum.
313			*/
314
315			if (!(r->crtype & REFRACTED)) { /* ray started in material */
316			VCOPY(v1, r->rdir);
317			n1[0] = -r->rdir[0]; n1[1] = -r->rdir[1]; n1[2] = -r->rdir[2];
318			} else {
319			/* find entry point */
320			for (entray = r; entray->rtype != REFRACTED;
321			entray = entray->parent)
322			;
323			entray = entray->parent;
324			if (entray->crtype & REFRACTED) /* too difficult */
325			return(0);
326			VCOPY(v1, entray->rdir);
327			VCOPY(n1, entray->ron);
328			}
329			VCOPY(v2, vt); /* exiting ray */
330			VCOPY(n2, r->ron);
331
332			/* first order dispersion approx. */
333			dtmp1 = DOT(n1, v1);
334			dtmp2 = DOT(n2, v2);
335			for (i = 0; i < 3; i++)
336			dv[i] = v1[i] + v2[i] - n1[i]/dtmp1 - n2[i]/dtmp2;
337
338			if (DOT(dv, dv) <= FTINY) /* null effect */
339			return(0);
340			/* compute plane normal */
341			fcross(v2Xdv, v2, dv);
342			v2Xdvv2Xdv = DOT(v2Xdv, v2Xdv);
343
344			/* check sources */
345	greg	1.7	initsrcindex(&si);
346			while (srcray(&sray, r, &si)) {
347	greg	1.1
348	greg	1.7	if (DOT(sray.rdir, v2) < MINCOS)
349	greg	1.1	continue; /* bad source */
350			/* adjust source ray */
351
352			dtmp1 = DOT(v2Xdv, sray.rdir) / v2Xdvv2Xdv;
353			sray.rdir[0] -= dtmp1 * v2Xdv[0];
354			sray.rdir[1] -= dtmp1 * v2Xdv[1];
355			sray.rdir[2] -= dtmp1 * v2Xdv[2];
356
357			l1 = lambda(m, v2, dv, sray.rdir); /* mean lambda */
358
359			if (l1 > MAXLAMBDA \|\| l1 < MINLAMBDA) /* not visible */
360			continue;
361			/* trace source ray */
362	greg	2.12	copycolor(sray.cext, cet);
363			copycolor(sray.albedo, abt);
364	greg	1.1	normalize(sray.rdir);
365			rayvalue(&sray);
366	greg	1.2	if (bright(sray.rcol) <= FTINY) /* missed it */
367	greg	1.1	continue;
368
369			/*
370			* Compute spectral sum over diameter of source.
371			* First find directions for rays going to opposite
372			* sides of source, then compute wavelengths for each.
373			*/
374
375			fcross(vtmp1, v2Xdv, sray.rdir);
376	greg	1.7	dtmp1 = sqrt(si.dom / v2Xdvv2Xdv / PI);
377	greg	1.1
378			/* compute first ray */
379			for (i = 0; i < 3; i++)
380			vtmp2[i] = sray.rdir[i] + dtmp1*vtmp1[i];
381
382			l1 = lambda(m, v2, dv, vtmp2); /* first lambda */
383			if (l1 < 0)
384			continue;
385			/* compute second ray */
386			for (i = 0; i < 3; i++)
387			vtmp2[i] = sray.rdir[i] - dtmp1*vtmp1[i];
388
389			l2 = lambda(m, v2, dv, vtmp2); /* second lambda */
390			if (l2 < 0)
391			continue;
392			/* compute color from spectrum */
393			if (l1 < l2)
394			spec_rgb(ctmp, l1, l2);
395			else
396			spec_rgb(ctmp, l2, l1);
397			multcolor(ctmp, sray.rcol);
398			scalecolor(ctmp, tr);
399			addcolor(r->rcol, ctmp);
400			success++;
401			}
402			return(success);
403			}
404
405
406			static int
407			lambda(m, v2, dv, lr) /* compute lambda for material */
408			register OBJREC *m;
409			FVECT v2, dv, lr;
410			{
411			FVECT lrXdv, v2Xlr;
412			double dtmp, denom;
413			int i;
414
415			fcross(lrXdv, lr, dv);
416			for (i = 0; i < 3; i++)
417			if (lrXdv[i] > FTINY \|\| lrXdv[i] < -FTINY)
418			break;
419			if (i >= 3)
420			return(-1);
421
422			fcross(v2Xlr, v2, lr);
423
424			dtmp = m->oargs.farg[4] / MLAMBDA;
425			denom = dtmp + v2Xlr[i]/lrXdv[i] * (m->oargs.farg[3] + dtmp);
426
427			if (denom < FTINY)
428			return(-1);
429
430			return(m->oargs.farg[4] / denom);
431			}
432
433			#endif /* DISPERSE */