ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/dielectric.c
Revision: 2.24
Committed: Sat May 10 17:43:01 2014 UTC (9 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad4R2, rad4R2P1
Changes since 2.23: +9 -8 lines
Log Message:
Fixed virtual distance so as not to undermine ambient calculation

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: dielectric.c,v 2.23 2013/08/07 05:10:09 greg 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 int
68 m_dielectric( /* color a ray which hit a dielectric interface */
69 OBJREC *m,
70 RAY *r
71 )
72 {
73 double cos1, cos2, nratio;
74 COLOR ctrans;
75 COLOR talb;
76 int hastexture;
77 double transdist=0, transtest=0;
78 double mirdist=0, mirtest=0;
79 int flatsurface;
80 double refl, trans;
81 FVECT dnorm;
82 double d1, d2;
83 RAY p;
84 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 = r->ro != NULL && isflat(r->ro->otype) &&
98 !hastexture | (r->crtype & AMBIENT);
99
100 /* index of refraction */
101 if (m->otype == MAT_DIELECTRIC)
102 nratio = m->oargs.farg[3] + m->oargs.farg[4]/MLAMBDA;
103 else
104 nratio = m->oargs.farg[3] / m->oargs.farg[7];
105
106 if (cos1 < 0.0) { /* inside */
107 hastexture = -hastexture;
108 cos1 = -cos1;
109 dnorm[0] = -dnorm[0];
110 dnorm[1] = -dnorm[1];
111 dnorm[2] = -dnorm[2];
112 setcolor(r->cext, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
113 -mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
114 -mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
115 setcolor(r->albedo, 0., 0., 0.);
116 r->gecc = 0.;
117 if (m->otype == MAT_INTERFACE) {
118 setcolor(ctrans,
119 -mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
120 -mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
121 -mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
122 setcolor(talb, 0., 0., 0.);
123 } else {
124 copycolor(ctrans, cextinction);
125 copycolor(talb, salbedo);
126 }
127 } else { /* outside */
128 nratio = 1.0 / nratio;
129
130 setcolor(ctrans, -mylog(m->oargs.farg[0]*colval(r->pcol,RED)),
131 -mylog(m->oargs.farg[1]*colval(r->pcol,GRN)),
132 -mylog(m->oargs.farg[2]*colval(r->pcol,BLU)));
133 setcolor(talb, 0., 0., 0.);
134 if (m->otype == MAT_INTERFACE) {
135 setcolor(r->cext,
136 -mylog(m->oargs.farg[4]*colval(r->pcol,RED)),
137 -mylog(m->oargs.farg[5]*colval(r->pcol,GRN)),
138 -mylog(m->oargs.farg[6]*colval(r->pcol,BLU)));
139 setcolor(r->albedo, 0., 0., 0.);
140 r->gecc = 0.;
141 }
142 }
143
144 d2 = 1.0 - nratio*nratio*(1.0 - cos1*cos1); /* compute cos theta2 */
145
146 if (d2 < FTINY) /* total reflection */
147
148 refl = 1.0;
149
150 else { /* refraction occurs */
151 /* compute Fresnel's equations */
152 cos2 = sqrt(d2);
153 d1 = cos1;
154 d2 = nratio*cos2;
155 d1 = (d1 - d2) / (d1 + d2);
156 refl = d1 * d1;
157
158 d1 = 1.0 / cos1;
159 d2 = nratio / cos2;
160 d1 = (d1 - d2) / (d1 + d2);
161 refl += d1 * d1;
162
163 refl *= 0.5;
164 trans = 1.0 - refl;
165
166 trans *= nratio*nratio; /* solid angle ratio */
167
168 setcolor(p.rcoef, trans, trans, trans);
169
170 if (rayorigin(&p, REFRACTED, r, p.rcoef) == 0) {
171
172 /* compute refracted ray */
173 d1 = nratio*cos1 - cos2;
174 for (i = 0; i < 3; i++)
175 p.rdir[i] = nratio*r->rdir[i] + d1*dnorm[i];
176 /* accidental reflection? */
177 if (hastexture &&
178 DOT(p.rdir,r->ron)*hastexture >= -FTINY) {
179 d1 *= (double)hastexture;
180 for (i = 0; i < 3; i++) /* ignore texture */
181 p.rdir[i] = nratio*r->rdir[i] +
182 d1*r->ron[i];
183 normalize(p.rdir); /* not exact */
184 } else
185 checknorm(p.rdir);
186 #ifdef DISPERSE
187 if (m->otype != MAT_DIELECTRIC
188 || r->rod > 0.0
189 || r->crtype & SHADOW
190 || !directvis
191 || m->oargs.farg[4] == 0.0
192 || !disperse(m, r, p.rdir,
193 trans, ctrans, talb))
194 #endif
195 {
196 copycolor(p.cext, ctrans);
197 copycolor(p.albedo, talb);
198 rayvalue(&p);
199 multcolor(p.rcol, p.rcoef);
200 addcolor(r->rcol, p.rcol);
201 /* virtual distance */
202 if (flatsurface ||
203 (1.-FTINY <= nratio) &
204 (nratio <= 1.+FTINY)) {
205 transtest = 2*bright(p.rcol);
206 transdist = r->rot + p.rt;
207 }
208 }
209 }
210 }
211 setcolor(p.rcoef, refl, refl, refl);
212
213 if (!(r->crtype & SHADOW) &&
214 rayorigin(&p, REFLECTED, r, p.rcoef) == 0) {
215
216 /* compute reflected ray */
217 VSUM(p.rdir, r->rdir, dnorm, 2.*cos1);
218 /* accidental penetration? */
219 if (hastexture && DOT(p.rdir,r->ron)*hastexture <= FTINY)
220 VSUM(p.rdir, r->rdir, r->ron, 2.*r->rod);
221 checknorm(p.rdir);
222 rayvalue(&p); /* reflected ray value */
223
224 multcolor(p.rcol, p.rcoef); /* color contribution */
225 addcolor(r->rcol, p.rcol);
226 /* virtual distance */
227 if (flatsurface) {
228 mirtest = 2*bright(p.rcol);
229 mirdist = r->rot + p.rt;
230 }
231 }
232 /* 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 /* rayvalue() computes absorption */
239 return(1);
240 }
241
242
243 #ifdef DISPERSE
244
245 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 {
255 RAY sray;
256 const RAY *entray;
257 FVECT v1, v2, n1, n2;
258 FVECT dv, v2Xdv;
259 double v2Xdvv2Xdv;
260 int success = 0;
261 SRCINDEX si;
262 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 = 1./DOT(n1, v1);
312 dtmp2 = 1./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 initsrcindex(&si);
324 while (srcray(&sray, r, &si)) {
325
326 if (DOT(sray.rdir, v2) < MINCOS)
327 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 copycolor(sray.cext, cet);
341 copycolor(sray.albedo, abt);
342 normalize(sray.rdir);
343 rayvalue(&sray);
344 if (bright(sray.rcol) <= FTINY) /* missed it */
345 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 dtmp1 = sqrt(si.dom / v2Xdvv2Xdv / PI);
355
356 /* compute first ray */
357 VSUM(vtmp2, sray.rdir, vtmp1, dtmp1);
358
359 l1 = lambda(m, v2, dv, vtmp2); /* first lambda */
360 if (l1 < 0)
361 continue;
362 /* compute second ray */
363 VSUM(vtmp2, sray.rdir, vtmp1, -dtmp1);
364
365 l2 = lambda(m, v2, dv, vtmp2); /* second lambda */
366 if (l2 < 0)
367 continue;
368 /* compute color from spectrum */
369 if (l1 < l2)
370 spec_rgb(ctmp, l1, l2);
371 else
372 spec_rgb(ctmp, l2, l1);
373 multcolor(ctmp, sray.rcol);
374 scalecolor(ctmp, tr);
375 addcolor(r->rcol, ctmp);
376 success++;
377 }
378 return(success);
379 }
380
381
382 static int
383 lambda( /* compute lambda for material */
384 OBJREC *m,
385 FVECT v2,
386 FVECT dv,
387 FVECT lr
388 )
389 {
390 FVECT lrXdv, v2Xlr;
391 double dtmp, denom;
392 int i;
393
394 fcross(lrXdv, lr, dv);
395 for (i = 0; i < 3; i++)
396 if ((lrXdv[i] > FTINY) | (lrXdv[i] < -FTINY))
397 break;
398 if (i >= 3)
399 return(-1);
400
401 fcross(v2Xlr, v2, lr);
402
403 dtmp = m->oargs.farg[4] / MLAMBDA;
404 denom = dtmp + v2Xlr[i]/lrXdv[i] * (m->oargs.farg[3] + dtmp);
405
406 if (denom < FTINY)
407 return(-1);
408
409 return(m->oargs.farg[4] / denom);
410 }
411
412 #endif /* DISPERSE */