ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/dielectric.c
Revision: 2.22
Committed: Mon Oct 25 22:57:45 2010 UTC (13 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R1
Changes since 2.21: +9 -13 lines
Log Message:
Changed to FVECT macros where appropriate

File Contents

# User Rev Content
1 greg 1.1 #ifndef lint
2 greg 2.22 static const char RCSid[] = "$Id: dielectric.c,v 2.21 2010/09/26 15:51:15 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 greg 2.21 } else
184     checknorm(p.rdir);
185 greg 1.1 #ifdef DISPERSE
186     if (m->otype != MAT_DIELECTRIC
187     || r->rod > 0.0
188     || r->crtype & SHADOW
189 greg 2.3 || !directvis
190 greg 1.1 || m->oargs.farg[4] == 0.0
191 greg 2.12 || !disperse(m, r, p.rdir,
192     trans, ctrans, talb))
193 greg 1.1 #endif
194     {
195 greg 2.9 copycolor(p.cext, ctrans);
196 greg 2.11 copycolor(p.albedo, talb);
197 greg 1.1 rayvalue(&p);
198 greg 2.20 multcolor(p.rcol, p.rcoef);
199 greg 1.1 addcolor(r->rcol, p.rcol);
200 greg 2.19 /* virtual distance */
201     if (flatsurface ||
202     (1.-FTINY <= nratio &&
203     nratio <= 1.+FTINY)) {
204     transtest = 2*bright(p.rcol);
205     transdist = r->rot + p.rt;
206     }
207 greg 1.1 }
208     }
209     }
210 greg 2.20 setcolor(p.rcoef, refl, refl, refl);
211    
212 greg 1.1 if (!(r->crtype & SHADOW) &&
213 greg 2.20 rayorigin(&p, REFLECTED, r, p.rcoef) == 0) {
214 greg 1.1
215     /* compute reflected ray */
216 greg 2.22 VSUM(p.rdir, r->rdir, dnorm, 2.*cos1);
217 gwlarson 2.14 /* accidental penetration? */
218     if (hastexture && DOT(p.rdir,r->ron)*hastexture <= FTINY)
219 greg 2.22 VSUM(p.rdir, r->rdir, r->ron, 2.*r->rod);
220 greg 2.21 checknorm(p.rdir);
221 greg 1.1 rayvalue(&p); /* reflected ray value */
222    
223 greg 2.20 multcolor(p.rcol, p.rcoef); /* color contribution */
224 greg 1.1 addcolor(r->rcol, p.rcol);
225 greg 2.19 /* virtual distance */
226     if (flatsurface) {
227     mirtest = 2*bright(p.rcol);
228     mirdist = r->rot + p.rt;
229     }
230 greg 1.1 }
231 greg 2.19 /* check distance to return */
232     d1 = bright(r->rcol);
233     if (transtest > d1)
234     r->rt = transdist;
235     else if (mirtest > d1)
236     r->rt = mirdist;
237 greg 2.9 /* rayvalue() computes absorption */
238 greg 2.7 return(1);
239 greg 1.1 }
240    
241    
242     #ifdef DISPERSE
243    
244 schorsch 2.18 static int
245     disperse( /* check light sources for dispersion */
246     OBJREC *m,
247     RAY *r,
248     FVECT vt,
249     double tr,
250     COLOR cet,
251     COLOR abt
252     )
253 greg 1.1 {
254 greg 2.20 RAY sray;
255     const RAY *entray;
256 greg 1.1 FVECT v1, v2, n1, n2;
257     FVECT dv, v2Xdv;
258     double v2Xdvv2Xdv;
259 greg 1.7 int success = 0;
260     SRCINDEX si;
261 greg 1.1 FVECT vtmp1, vtmp2;
262     double dtmp1, dtmp2;
263     int l1, l2;
264     COLOR ctmp;
265     int i;
266    
267     /*
268     * This routine computes dispersion to the first order using
269     * the following assumptions:
270     *
271     * 1) The dependency of the index of refraction on wavelength
272     * is approximated by Hartmann's equation with lambda0
273     * equal to zero.
274     * 2) The entry and exit locations are constant with respect
275     * to dispersion.
276     *
277     * The second assumption permits us to model dispersion without
278     * having to sample refracted directions. We assume that the
279     * geometry inside the material is constant, and concern ourselves
280     * only with the relationship between the entering and exiting ray.
281     * We compute the first derivatives of the entering and exiting
282     * refraction with respect to the index of refraction. This
283     * is then used in a first order Taylor series to determine the
284     * index of refraction necessary to send the exiting ray to each
285     * light source.
286     * If an exiting ray hits a light source within the refraction
287     * boundaries, we sum all the frequencies over the disc of the
288     * light source to determine the resulting color. A smaller light
289     * source will therefore exhibit a sharper spectrum.
290     */
291    
292     if (!(r->crtype & REFRACTED)) { /* ray started in material */
293     VCOPY(v1, r->rdir);
294     n1[0] = -r->rdir[0]; n1[1] = -r->rdir[1]; n1[2] = -r->rdir[2];
295     } else {
296     /* find entry point */
297     for (entray = r; entray->rtype != REFRACTED;
298     entray = entray->parent)
299     ;
300     entray = entray->parent;
301     if (entray->crtype & REFRACTED) /* too difficult */
302     return(0);
303     VCOPY(v1, entray->rdir);
304     VCOPY(n1, entray->ron);
305     }
306     VCOPY(v2, vt); /* exiting ray */
307     VCOPY(n2, r->ron);
308    
309     /* first order dispersion approx. */
310 greg 2.22 dtmp1 = 1./DOT(n1, v1);
311     dtmp2 = 1./DOT(n2, v2);
312 greg 1.1 for (i = 0; i < 3; i++)
313 greg 2.22 dv[i] = v1[i] + v2[i] - n1[i]*dtmp1 - n2[i]*dtmp2;
314 greg 1.1
315     if (DOT(dv, dv) <= FTINY) /* null effect */
316     return(0);
317     /* compute plane normal */
318     fcross(v2Xdv, v2, dv);
319     v2Xdvv2Xdv = DOT(v2Xdv, v2Xdv);
320    
321     /* check sources */
322 greg 1.7 initsrcindex(&si);
323     while (srcray(&sray, r, &si)) {
324 greg 1.1
325 greg 1.7 if (DOT(sray.rdir, v2) < MINCOS)
326 greg 1.1 continue; /* bad source */
327     /* adjust source ray */
328    
329     dtmp1 = DOT(v2Xdv, sray.rdir) / v2Xdvv2Xdv;
330     sray.rdir[0] -= dtmp1 * v2Xdv[0];
331     sray.rdir[1] -= dtmp1 * v2Xdv[1];
332     sray.rdir[2] -= dtmp1 * v2Xdv[2];
333    
334     l1 = lambda(m, v2, dv, sray.rdir); /* mean lambda */
335    
336     if (l1 > MAXLAMBDA || l1 < MINLAMBDA) /* not visible */
337     continue;
338     /* trace source ray */
339 greg 2.12 copycolor(sray.cext, cet);
340     copycolor(sray.albedo, abt);
341 greg 1.1 normalize(sray.rdir);
342     rayvalue(&sray);
343 greg 1.2 if (bright(sray.rcol) <= FTINY) /* missed it */
344 greg 1.1 continue;
345    
346     /*
347     * Compute spectral sum over diameter of source.
348     * First find directions for rays going to opposite
349     * sides of source, then compute wavelengths for each.
350     */
351    
352     fcross(vtmp1, v2Xdv, sray.rdir);
353 greg 1.7 dtmp1 = sqrt(si.dom / v2Xdvv2Xdv / PI);
354 greg 1.1
355     /* compute first ray */
356 greg 2.22 VSUM(vtmp2, sray.rdir, vtmp1, dtmp1);
357 greg 1.1
358     l1 = lambda(m, v2, dv, vtmp2); /* first lambda */
359     if (l1 < 0)
360     continue;
361     /* compute second ray */
362 greg 2.22 VSUM(vtmp2, sray.rdir, vtmp1, -dtmp1);
363 greg 1.1
364     l2 = lambda(m, v2, dv, vtmp2); /* second lambda */
365     if (l2 < 0)
366     continue;
367     /* compute color from spectrum */
368     if (l1 < l2)
369     spec_rgb(ctmp, l1, l2);
370     else
371     spec_rgb(ctmp, l2, l1);
372     multcolor(ctmp, sray.rcol);
373     scalecolor(ctmp, tr);
374     addcolor(r->rcol, ctmp);
375     success++;
376     }
377     return(success);
378     }
379    
380    
381     static int
382 schorsch 2.18 lambda( /* compute lambda for material */
383     register OBJREC *m,
384     FVECT v2,
385     FVECT dv,
386     FVECT lr
387     )
388 greg 1.1 {
389     FVECT lrXdv, v2Xlr;
390     double dtmp, denom;
391     int i;
392    
393     fcross(lrXdv, lr, dv);
394     for (i = 0; i < 3; i++)
395 greg 2.22 if ((lrXdv[i] > FTINY) | (lrXdv[i] < -FTINY))
396 greg 1.1 break;
397     if (i >= 3)
398     return(-1);
399    
400     fcross(v2Xlr, v2, lr);
401    
402     dtmp = m->oargs.farg[4] / MLAMBDA;
403     denom = dtmp + v2Xlr[i]/lrXdv[i] * (m->oargs.farg[3] + dtmp);
404    
405     if (denom < FTINY)
406     return(-1);
407    
408     return(m->oargs.farg[4] / denom);
409     }
410    
411     #endif /* DISPERSE */