ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/dielectric.c
Revision: 2.19
Committed: Thu Sep 9 06:46:07 2004 UTC (19 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R6, rad3R6P1
Changes since 2.18: +24 -3 lines
Log Message:
Corrected effective ray length calculation

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: dielectric.c,v 2.18 2004/03/30 16:13:01 schorsch Exp $";
3 #endif
4 /*
5 * dielectric.c - shading function for transparent materials.
6 */
7
8 #include "copyright.h"
9
10 #include "ray.h"
11 #include "otypes.h"
12 #include "rtotypes.h"
13
14 #ifdef DISPERSE
15 #include "source.h"
16 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 #endif
19
20 static double mylog(double x);
21
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
54 static double
55 mylog( /* special log for extinction coefficients */
56 double x
57 )
58 {
59 if (x < 1e-40)
60 return(-100.);
61 if (x >= 1.)
62 return(0.);
63 return(log(x));
64 }
65
66
67 extern int
68 m_dielectric( /* color a ray which hit a dielectric interface */
69 OBJREC *m,
70 register RAY *r
71 )
72 {
73 double cos1, cos2, nratio;
74 COLOR ctrans;
75 COLOR talb;
76 int hastexture;
77 double transdist, transtest=0;
78 double mirdist, mirtest=0;
79 int flatsurface;
80 double refl, trans;
81 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 if ( (hastexture = DOT(r->pert,r->pert) > FTINY*FTINY) )
92 cos1 = raynormal(dnorm, r); /* perturb normal */
93 else {
94 VCOPY(dnorm, r->ron);
95 cos1 = r->rod;
96 }
97 flatsurface = !hastexture && r->ro != NULL && isflat(r->ro->otype);
98
99 /* 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 hastexture = -hastexture;
107 cos1 = -cos1;
108 dnorm[0] = -dnorm[0];
109 dnorm[1] = -dnorm[1];
110 dnorm[2] = -dnorm[2];
111 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 setcolor(r->albedo, 0., 0., 0.);
115 r->gecc = 0.;
116 if (m->otype == MAT_INTERFACE) {
117 setcolor(ctrans,
118 -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 setcolor(talb, 0., 0., 0.);
122 } else {
123 copycolor(ctrans, cextinction);
124 copycolor(talb, salbedo);
125 }
126 } else { /* outside */
127 nratio = 1.0 / nratio;
128
129 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 setcolor(talb, 0., 0., 0.);
133 if (m->otype == MAT_INTERFACE) {
134 setcolor(r->cext,
135 -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 setcolor(r->albedo, 0., 0., 0.);
139 r->gecc = 0.;
140 }
141 }
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 refl *= 0.5;
163 trans = 1.0 - refl;
164
165 trans *= nratio*nratio; /* solid angle ratio */
166
167 if (rayorigin(&p, r, REFRACTED, trans) == 0) {
168
169 /* compute refracted ray */
170 d1 = nratio*cos1 - cos2;
171 for (i = 0; i < 3; i++)
172 p.rdir[i] = nratio*r->rdir[i] + d1*dnorm[i];
173 /* accidental reflection? */
174 if (hastexture &&
175 DOT(p.rdir,r->ron)*hastexture >= -FTINY) {
176 d1 *= (double)hastexture;
177 for (i = 0; i < 3; i++) /* ignore texture */
178 p.rdir[i] = nratio*r->rdir[i] +
179 d1*r->ron[i];
180 normalize(p.rdir); /* not exact */
181 }
182 #ifdef DISPERSE
183 if (m->otype != MAT_DIELECTRIC
184 || r->rod > 0.0
185 || r->crtype & SHADOW
186 || !directvis
187 || m->oargs.farg[4] == 0.0
188 || !disperse(m, r, p.rdir,
189 trans, ctrans, talb))
190 #endif
191 {
192 copycolor(p.cext, ctrans);
193 copycolor(p.albedo, talb);
194 rayvalue(&p);
195 scalecolor(p.rcol, trans);
196 addcolor(r->rcol, p.rcol);
197 /* virtual distance */
198 if (flatsurface ||
199 (1.-FTINY <= nratio &&
200 nratio <= 1.+FTINY)) {
201 transtest = 2*bright(p.rcol);
202 transdist = r->rot + p.rt;
203 }
204 }
205 }
206 }
207
208 if (!(r->crtype & SHADOW) &&
209 rayorigin(&p, r, REFLECTED, refl) == 0) {
210
211 /* compute reflected ray */
212 for (i = 0; i < 3; i++)
213 p.rdir[i] = r->rdir[i] + 2.0*cos1*dnorm[i];
214 /* accidental penetration? */
215 if (hastexture && DOT(p.rdir,r->ron)*hastexture <= FTINY)
216 for (i = 0; i < 3; i++) /* ignore texture */
217 p.rdir[i] = r->rdir[i] + 2.0*r->rod*r->ron[i];
218
219 rayvalue(&p); /* reflected ray value */
220
221 scalecolor(p.rcol, refl); /* color contribution */
222 addcolor(r->rcol, p.rcol);
223 /* virtual distance */
224 if (flatsurface) {
225 mirtest = 2*bright(p.rcol);
226 mirdist = r->rot + p.rt;
227 }
228 }
229 /* check distance to return */
230 d1 = bright(r->rcol);
231 if (transtest > d1)
232 r->rt = transdist;
233 else if (mirtest > d1)
234 r->rt = mirdist;
235 /* rayvalue() computes absorption */
236 return(1);
237 }
238
239
240 #ifdef DISPERSE
241
242 static int
243 disperse( /* check light sources for dispersion */
244 OBJREC *m,
245 RAY *r,
246 FVECT vt,
247 double tr,
248 COLOR cet,
249 COLOR abt
250 )
251 {
252 RAY sray, *entray;
253 FVECT v1, v2, n1, n2;
254 FVECT dv, v2Xdv;
255 double v2Xdvv2Xdv;
256 int success = 0;
257 SRCINDEX si;
258 FVECT vtmp1, vtmp2;
259 double dtmp1, dtmp2;
260 int l1, l2;
261 COLOR ctmp;
262 int i;
263
264 /*
265 * This routine computes dispersion to the first order using
266 * the following assumptions:
267 *
268 * 1) The dependency of the index of refraction on wavelength
269 * is approximated by Hartmann's equation with lambda0
270 * equal to zero.
271 * 2) The entry and exit locations are constant with respect
272 * to dispersion.
273 *
274 * The second assumption permits us to model dispersion without
275 * having to sample refracted directions. We assume that the
276 * geometry inside the material is constant, and concern ourselves
277 * only with the relationship between the entering and exiting ray.
278 * We compute the first derivatives of the entering and exiting
279 * refraction with respect to the index of refraction. This
280 * is then used in a first order Taylor series to determine the
281 * index of refraction necessary to send the exiting ray to each
282 * light source.
283 * If an exiting ray hits a light source within the refraction
284 * boundaries, we sum all the frequencies over the disc of the
285 * light source to determine the resulting color. A smaller light
286 * source will therefore exhibit a sharper spectrum.
287 */
288
289 if (!(r->crtype & REFRACTED)) { /* ray started in material */
290 VCOPY(v1, r->rdir);
291 n1[0] = -r->rdir[0]; n1[1] = -r->rdir[1]; n1[2] = -r->rdir[2];
292 } else {
293 /* find entry point */
294 for (entray = r; entray->rtype != REFRACTED;
295 entray = entray->parent)
296 ;
297 entray = entray->parent;
298 if (entray->crtype & REFRACTED) /* too difficult */
299 return(0);
300 VCOPY(v1, entray->rdir);
301 VCOPY(n1, entray->ron);
302 }
303 VCOPY(v2, vt); /* exiting ray */
304 VCOPY(n2, r->ron);
305
306 /* first order dispersion approx. */
307 dtmp1 = DOT(n1, v1);
308 dtmp2 = DOT(n2, v2);
309 for (i = 0; i < 3; i++)
310 dv[i] = v1[i] + v2[i] - n1[i]/dtmp1 - n2[i]/dtmp2;
311
312 if (DOT(dv, dv) <= FTINY) /* null effect */
313 return(0);
314 /* compute plane normal */
315 fcross(v2Xdv, v2, dv);
316 v2Xdvv2Xdv = DOT(v2Xdv, v2Xdv);
317
318 /* check sources */
319 initsrcindex(&si);
320 while (srcray(&sray, r, &si)) {
321
322 if (DOT(sray.rdir, v2) < MINCOS)
323 continue; /* bad source */
324 /* adjust source ray */
325
326 dtmp1 = DOT(v2Xdv, sray.rdir) / v2Xdvv2Xdv;
327 sray.rdir[0] -= dtmp1 * v2Xdv[0];
328 sray.rdir[1] -= dtmp1 * v2Xdv[1];
329 sray.rdir[2] -= dtmp1 * v2Xdv[2];
330
331 l1 = lambda(m, v2, dv, sray.rdir); /* mean lambda */
332
333 if (l1 > MAXLAMBDA || l1 < MINLAMBDA) /* not visible */
334 continue;
335 /* trace source ray */
336 copycolor(sray.cext, cet);
337 copycolor(sray.albedo, abt);
338 normalize(sray.rdir);
339 rayvalue(&sray);
340 if (bright(sray.rcol) <= FTINY) /* missed it */
341 continue;
342
343 /*
344 * Compute spectral sum over diameter of source.
345 * First find directions for rays going to opposite
346 * sides of source, then compute wavelengths for each.
347 */
348
349 fcross(vtmp1, v2Xdv, sray.rdir);
350 dtmp1 = sqrt(si.dom / v2Xdvv2Xdv / PI);
351
352 /* compute first ray */
353 for (i = 0; i < 3; i++)
354 vtmp2[i] = sray.rdir[i] + dtmp1*vtmp1[i];
355
356 l1 = lambda(m, v2, dv, vtmp2); /* first lambda */
357 if (l1 < 0)
358 continue;
359 /* compute second ray */
360 for (i = 0; i < 3; i++)
361 vtmp2[i] = sray.rdir[i] - dtmp1*vtmp1[i];
362
363 l2 = lambda(m, v2, dv, vtmp2); /* second lambda */
364 if (l2 < 0)
365 continue;
366 /* compute color from spectrum */
367 if (l1 < l2)
368 spec_rgb(ctmp, l1, l2);
369 else
370 spec_rgb(ctmp, l2, l1);
371 multcolor(ctmp, sray.rcol);
372 scalecolor(ctmp, tr);
373 addcolor(r->rcol, ctmp);
374 success++;
375 }
376 return(success);
377 }
378
379
380 static int
381 lambda( /* compute lambda for material */
382 register OBJREC *m,
383 FVECT v2,
384 FVECT dv,
385 FVECT lr
386 )
387 {
388 FVECT lrXdv, v2Xlr;
389 double dtmp, denom;
390 int i;
391
392 fcross(lrXdv, lr, dv);
393 for (i = 0; i < 3; i++)
394 if (lrXdv[i] > FTINY || lrXdv[i] < -FTINY)
395 break;
396 if (i >= 3)
397 return(-1);
398
399 fcross(v2Xlr, v2, lr);
400
401 dtmp = m->oargs.farg[4] / MLAMBDA;
402 denom = dtmp + v2Xlr[i]/lrXdv[i] * (m->oargs.farg[3] + dtmp);
403
404 if (denom < FTINY)
405 return(-1);
406
407 return(m->oargs.farg[4] / denom);
408 }
409
410 #endif /* DISPERSE */