ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/dielectric.c
Revision: 2.20
Committed: Tue Apr 19 01:15:06 2005 UTC (19 years ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad4R0, rad3R8, rad3R9
Changes since 2.19: +11 -7 lines
Log Message:
Extensive changes to enable rtrace -oTW option for tracking ray contributions

File Contents

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