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, 7 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

# User Rev Content
1 greg 1.1 #ifndef lint
2 greg 2.19 static const char RCSid[] = "$Id: dielectric.c,v 2.18 2004/03/30 16:13:01 schorsch 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 gwlarson 2.13 if (rayorigin(&p, r, REFRACTED, trans) == 0) {
168 greg 1.1
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 gwlarson 2.14 /* 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 greg 1.1 #ifdef DISPERSE
183     if (m->otype != MAT_DIELECTRIC
184     || r->rod > 0.0
185     || r->crtype & SHADOW
186 greg 2.3 || !directvis
187 greg 1.1 || m->oargs.farg[4] == 0.0
188 greg 2.12 || !disperse(m, r, p.rdir,
189     trans, ctrans, talb))
190 greg 1.1 #endif
191     {
192 greg 2.9 copycolor(p.cext, ctrans);
193 greg 2.11 copycolor(p.albedo, talb);
194 greg 1.1 rayvalue(&p);
195     scalecolor(p.rcol, trans);
196     addcolor(r->rcol, p.rcol);
197 greg 2.19 /* 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 greg 1.1 }
205     }
206     }
207    
208     if (!(r->crtype & SHADOW) &&
209 gwlarson 2.13 rayorigin(&p, r, REFLECTED, refl) == 0) {
210 greg 1.1
211     /* compute reflected ray */
212     for (i = 0; i < 3; i++)
213     p.rdir[i] = r->rdir[i] + 2.0*cos1*dnorm[i];
214 gwlarson 2.14 /* 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 greg 1.1
219     rayvalue(&p); /* reflected ray value */
220    
221     scalecolor(p.rcol, refl); /* color contribution */
222     addcolor(r->rcol, p.rcol);
223 greg 2.19 /* virtual distance */
224     if (flatsurface) {
225     mirtest = 2*bright(p.rcol);
226     mirdist = r->rot + p.rt;
227     }
228 greg 1.1 }
229 greg 2.19 /* 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 greg 2.9 /* rayvalue() computes absorption */
236 greg 2.7 return(1);
237 greg 1.1 }
238    
239    
240     #ifdef DISPERSE
241    
242 schorsch 2.18 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 greg 1.1 {
252     RAY sray, *entray;
253     FVECT v1, v2, n1, n2;
254     FVECT dv, v2Xdv;
255     double v2Xdvv2Xdv;
256 greg 1.7 int success = 0;
257     SRCINDEX si;
258 greg 1.1 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 greg 1.7 initsrcindex(&si);
320     while (srcray(&sray, r, &si)) {
321 greg 1.1
322 greg 1.7 if (DOT(sray.rdir, v2) < MINCOS)
323 greg 1.1 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 greg 2.12 copycolor(sray.cext, cet);
337     copycolor(sray.albedo, abt);
338 greg 1.1 normalize(sray.rdir);
339     rayvalue(&sray);
340 greg 1.2 if (bright(sray.rcol) <= FTINY) /* missed it */
341 greg 1.1 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 greg 1.7 dtmp1 = sqrt(si.dom / v2Xdvv2Xdv / PI);
351 greg 1.1
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 schorsch 2.18 lambda( /* compute lambda for material */
382     register OBJREC *m,
383     FVECT v2,
384     FVECT dv,
385     FVECT lr
386     )
387 greg 1.1 {
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 */