ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/dielectric.c
Revision: 2.12
Committed: Thu Jul 11 15:14:05 1996 UTC (27 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.11: +14 -3 lines
Log Message:
fixed bug in computation of absorption and dispersion routines

File Contents

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