ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/dielectric.c
Revision: 2.13
Committed: Wed Jun 17 12:53:07 1998 UTC (25 years, 10 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 2.12: +4 -13 lines
Log Message:
improved the way ray weights are computed for participating media

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