ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/dielectric.c
Revision: 2.14
Committed: Fri Jun 19 12:01:15 1998 UTC (25 years, 10 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 2.13: +21 -2 lines
Log Message:
added tests for accidental reflection/penetration

File Contents

# User Rev Content
1 gwlarson 2.13 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2 greg 1.1
3     #ifndef lint
4 gwlarson 2.13 static char SCCSid[] = "$SunId$ SGI";
5 greg 1.1 #endif
6    
7     /*
8     * dielectric.c - shading function for transparent materials.
9     *
10     * 9/6/85
11     */
12    
13     #include "ray.h"
14    
15     #include "otypes.h"
16    
17     #ifdef DISPERSE
18     #include "source.h"
19 greg 2.5 static disperse();
20 greg 2.6 static int lambda();
21 greg 1.1 #endif
22    
23     /*
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 extern COLOR cextinction; /* global coefficient of extinction */
54 greg 2.11 extern COLOR salbedo; /* global scattering albedo */
55 greg 1.1
56 greg 2.9
57 greg 2.10 static double
58     mylog(x) /* special log for extinction coefficients */
59     double x;
60     {
61     if (x < 1e-40)
62     return(-100.);
63     if (x >= 1.)
64     return(0.);
65     return(log(x));
66     }
67    
68    
69 greg 2.9 m_dielectric(m, r) /* color a ray which hit a dielectric interface */
70 greg 1.1 OBJREC *m;
71     register RAY *r;
72     {
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 1.5 double refl, trans;
78 greg 1.1 FVECT dnorm;
79     double d1, d2;
80     RAY p;
81     register int i;
82    
83     if (m->oargs.nfargs != (m->otype==MAT_DIELECTRIC ? 5 : 8))
84     objerror(m, USER, "bad arguments");
85    
86     raytexture(r, m->omod); /* get modifiers */
87    
88 gwlarson 2.14 if (hastexture = DOT(r->pert,r->pert) > FTINY*FTINY)
89     cos1 = raynormal(dnorm, r); /* perturb normal */
90     else {
91     VCOPY(dnorm, r->ron);
92     cos1 = r->rod;
93     }
94 greg 1.1 /* index of refraction */
95     if (m->otype == MAT_DIELECTRIC)
96     nratio = m->oargs.farg[3] + m->oargs.farg[4]/MLAMBDA;
97     else
98     nratio = m->oargs.farg[3] / m->oargs.farg[7];
99    
100     if (cos1 < 0.0) { /* inside */
101 gwlarson 2.14 hastexture = -hastexture;
102 greg 1.1 cos1 = -cos1;
103     dnorm[0] = -dnorm[0];
104     dnorm[1] = -dnorm[1];
105     dnorm[2] = -dnorm[2];
106 greg 2.10 setcolor(r->cext, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
107     -mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
108     -mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
109 greg 2.11 setcolor(r->albedo, 0., 0., 0.);
110 greg 2.9 r->gecc = 0.;
111     if (m->otype == MAT_INTERFACE) {
112     setcolor(ctrans,
113 greg 2.10 -mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
114     -mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
115     -mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
116 greg 2.11 setcolor(talb, 0., 0., 0.);
117 greg 2.9 } else {
118     copycolor(ctrans, cextinction);
119 greg 2.11 copycolor(talb, salbedo);
120 greg 2.9 }
121 greg 1.1 } else { /* outside */
122     nratio = 1.0 / nratio;
123 greg 2.9
124 greg 2.10 setcolor(ctrans, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
125     -mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
126     -mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
127 greg 2.11 setcolor(talb, 0., 0., 0.);
128 greg 2.9 if (m->otype == MAT_INTERFACE) {
129     setcolor(r->cext,
130 greg 2.10 -mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
131     -mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
132     -mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
133 greg 2.11 setcolor(r->albedo, 0., 0., 0.);
134 greg 2.9 r->gecc = 0.;
135     }
136 greg 1.1 }
137    
138     d2 = 1.0 - nratio*nratio*(1.0 - cos1*cos1); /* compute cos theta2 */
139    
140     if (d2 < FTINY) /* total reflection */
141    
142     refl = 1.0;
143    
144     else { /* refraction occurs */
145     /* compute Fresnel's equations */
146     cos2 = sqrt(d2);
147     d1 = cos1;
148     d2 = nratio*cos2;
149     d1 = (d1 - d2) / (d1 + d2);
150     refl = d1 * d1;
151    
152     d1 = 1.0 / cos1;
153     d2 = nratio / cos2;
154     d1 = (d1 - d2) / (d1 + d2);
155     refl += d1 * d1;
156    
157 greg 2.9 refl *= 0.5;
158 greg 1.1 trans = 1.0 - refl;
159    
160 gwlarson 2.13 if (rayorigin(&p, r, REFRACTED, trans) == 0) {
161 greg 1.1
162     /* compute refracted ray */
163     d1 = nratio*cos1 - cos2;
164     for (i = 0; i < 3; i++)
165     p.rdir[i] = nratio*r->rdir[i] + d1*dnorm[i];
166 gwlarson 2.14 /* accidental reflection? */
167     if (hastexture &&
168     DOT(p.rdir,r->ron)*hastexture >= -FTINY) {
169     d1 *= (double)hastexture;
170     for (i = 0; i < 3; i++) /* ignore texture */
171     p.rdir[i] = nratio*r->rdir[i] +
172     d1*r->ron[i];
173     normalize(p.rdir); /* not exact */
174     }
175 greg 1.1 #ifdef DISPERSE
176     if (m->otype != MAT_DIELECTRIC
177     || r->rod > 0.0
178     || r->crtype & SHADOW
179 greg 2.3 || !directvis
180 greg 1.1 || m->oargs.farg[4] == 0.0
181 greg 2.12 || !disperse(m, r, p.rdir,
182     trans, ctrans, talb))
183 greg 1.1 #endif
184     {
185 greg 2.9 copycolor(p.cext, ctrans);
186 greg 2.11 copycolor(p.albedo, talb);
187 greg 1.1 rayvalue(&p);
188     scalecolor(p.rcol, trans);
189     addcolor(r->rcol, p.rcol);
190 greg 2.4 if (nratio >= 1.0-FTINY && nratio <= 1.0+FTINY)
191     r->rt = r->rot + p.rt;
192 greg 1.1 }
193     }
194     }
195    
196     if (!(r->crtype & SHADOW) &&
197 gwlarson 2.13 rayorigin(&p, r, REFLECTED, refl) == 0) {
198 greg 1.1
199     /* compute reflected ray */
200     for (i = 0; i < 3; i++)
201     p.rdir[i] = r->rdir[i] + 2.0*cos1*dnorm[i];
202 gwlarson 2.14 /* accidental penetration? */
203     if (hastexture && DOT(p.rdir,r->ron)*hastexture <= FTINY)
204     for (i = 0; i < 3; i++) /* ignore texture */
205     p.rdir[i] = r->rdir[i] + 2.0*r->rod*r->ron[i];
206 greg 1.1
207     rayvalue(&p); /* reflected ray value */
208    
209     scalecolor(p.rcol, refl); /* color contribution */
210     addcolor(r->rcol, p.rcol);
211     }
212 greg 2.9 /* rayvalue() computes absorption */
213 greg 2.7 return(1);
214 greg 1.1 }
215    
216    
217     #ifdef DISPERSE
218    
219     static
220 greg 2.12 disperse(m, r, vt, tr, cet, abt) /* check light sources for dispersion */
221 greg 1.1 OBJREC *m;
222     RAY *r;
223     FVECT vt;
224     double tr;
225 greg 2.12 COLOR cet, abt;
226 greg 1.1 {
227     RAY sray, *entray;
228     FVECT v1, v2, n1, n2;
229     FVECT dv, v2Xdv;
230     double v2Xdvv2Xdv;
231 greg 1.7 int success = 0;
232     SRCINDEX si;
233 greg 1.1 FVECT vtmp1, vtmp2;
234     double dtmp1, dtmp2;
235     int l1, l2;
236     COLOR ctmp;
237     int i;
238    
239     /*
240     * This routine computes dispersion to the first order using
241     * the following assumptions:
242     *
243     * 1) The dependency of the index of refraction on wavelength
244     * is approximated by Hartmann's equation with lambda0
245     * equal to zero.
246     * 2) The entry and exit locations are constant with respect
247     * to dispersion.
248     *
249     * The second assumption permits us to model dispersion without
250     * having to sample refracted directions. We assume that the
251     * geometry inside the material is constant, and concern ourselves
252     * only with the relationship between the entering and exiting ray.
253     * We compute the first derivatives of the entering and exiting
254     * refraction with respect to the index of refraction. This
255     * is then used in a first order Taylor series to determine the
256     * index of refraction necessary to send the exiting ray to each
257     * light source.
258     * If an exiting ray hits a light source within the refraction
259     * boundaries, we sum all the frequencies over the disc of the
260     * light source to determine the resulting color. A smaller light
261     * source will therefore exhibit a sharper spectrum.
262     */
263    
264     if (!(r->crtype & REFRACTED)) { /* ray started in material */
265     VCOPY(v1, r->rdir);
266     n1[0] = -r->rdir[0]; n1[1] = -r->rdir[1]; n1[2] = -r->rdir[2];
267     } else {
268     /* find entry point */
269     for (entray = r; entray->rtype != REFRACTED;
270     entray = entray->parent)
271     ;
272     entray = entray->parent;
273     if (entray->crtype & REFRACTED) /* too difficult */
274     return(0);
275     VCOPY(v1, entray->rdir);
276     VCOPY(n1, entray->ron);
277     }
278     VCOPY(v2, vt); /* exiting ray */
279     VCOPY(n2, r->ron);
280    
281     /* first order dispersion approx. */
282     dtmp1 = DOT(n1, v1);
283     dtmp2 = DOT(n2, v2);
284     for (i = 0; i < 3; i++)
285     dv[i] = v1[i] + v2[i] - n1[i]/dtmp1 - n2[i]/dtmp2;
286    
287     if (DOT(dv, dv) <= FTINY) /* null effect */
288     return(0);
289     /* compute plane normal */
290     fcross(v2Xdv, v2, dv);
291     v2Xdvv2Xdv = DOT(v2Xdv, v2Xdv);
292    
293     /* check sources */
294 greg 1.7 initsrcindex(&si);
295     while (srcray(&sray, r, &si)) {
296 greg 1.1
297 greg 1.7 if (DOT(sray.rdir, v2) < MINCOS)
298 greg 1.1 continue; /* bad source */
299     /* adjust source ray */
300    
301     dtmp1 = DOT(v2Xdv, sray.rdir) / v2Xdvv2Xdv;
302     sray.rdir[0] -= dtmp1 * v2Xdv[0];
303     sray.rdir[1] -= dtmp1 * v2Xdv[1];
304     sray.rdir[2] -= dtmp1 * v2Xdv[2];
305    
306     l1 = lambda(m, v2, dv, sray.rdir); /* mean lambda */
307    
308     if (l1 > MAXLAMBDA || l1 < MINLAMBDA) /* not visible */
309     continue;
310     /* trace source ray */
311 greg 2.12 copycolor(sray.cext, cet);
312     copycolor(sray.albedo, abt);
313 greg 1.1 normalize(sray.rdir);
314     rayvalue(&sray);
315 greg 1.2 if (bright(sray.rcol) <= FTINY) /* missed it */
316 greg 1.1 continue;
317    
318     /*
319     * Compute spectral sum over diameter of source.
320     * First find directions for rays going to opposite
321     * sides of source, then compute wavelengths for each.
322     */
323    
324     fcross(vtmp1, v2Xdv, sray.rdir);
325 greg 1.7 dtmp1 = sqrt(si.dom / v2Xdvv2Xdv / PI);
326 greg 1.1
327     /* compute first ray */
328     for (i = 0; i < 3; i++)
329     vtmp2[i] = sray.rdir[i] + dtmp1*vtmp1[i];
330    
331     l1 = lambda(m, v2, dv, vtmp2); /* first lambda */
332     if (l1 < 0)
333     continue;
334     /* compute second ray */
335     for (i = 0; i < 3; i++)
336     vtmp2[i] = sray.rdir[i] - dtmp1*vtmp1[i];
337    
338     l2 = lambda(m, v2, dv, vtmp2); /* second lambda */
339     if (l2 < 0)
340     continue;
341     /* compute color from spectrum */
342     if (l1 < l2)
343     spec_rgb(ctmp, l1, l2);
344     else
345     spec_rgb(ctmp, l2, l1);
346     multcolor(ctmp, sray.rcol);
347     scalecolor(ctmp, tr);
348     addcolor(r->rcol, ctmp);
349     success++;
350     }
351     return(success);
352     }
353    
354    
355     static int
356     lambda(m, v2, dv, lr) /* compute lambda for material */
357     register OBJREC *m;
358     FVECT v2, dv, lr;
359     {
360     FVECT lrXdv, v2Xlr;
361     double dtmp, denom;
362     int i;
363    
364     fcross(lrXdv, lr, dv);
365     for (i = 0; i < 3; i++)
366     if (lrXdv[i] > FTINY || lrXdv[i] < -FTINY)
367     break;
368     if (i >= 3)
369     return(-1);
370    
371     fcross(v2Xlr, v2, lr);
372    
373     dtmp = m->oargs.farg[4] / MLAMBDA;
374     denom = dtmp + v2Xlr[i]/lrXdv[i] * (m->oargs.farg[3] + dtmp);
375    
376     if (denom < FTINY)
377     return(-1);
378    
379     return(m->oargs.farg[4] / denom);
380     }
381    
382     #endif /* DISPERSE */