home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 8 / CDASC08.ISO / NEWS / RADIANCE / SRC / RT / DIELECTR.C < prev    next >
C/C++ Source or Header  |  1993-10-07  |  9KB  |  324 lines
  1. /* Copyright (c) 1986 Regents of the University of California */
  2.  
  3. #ifndef lint
  4. static char SCCSid[] = "@(#)dielectric.c 2.2 10/2/92 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. #endif
  20.  
  21. /*
  22.  *     Explicit calculations for Fresnel's equation are performed,
  23.  *  but only one square root computation is necessary.
  24.  *     The index of refraction is given as a Hartmann equation
  25.  *  with lambda0 equal to zero.  If the slope of Hartmann's
  26.  *  equation is non-zero, the material disperses light upon
  27.  *  refraction.  This condition is examined on rays traced to
  28.  *  light sources.  If a ray is exiting a dielectric material, we
  29.  *  check the sources to see if any would cause bright color to be
  30.  *  directed to the viewer due to dispersion.  This gives colorful
  31.  *  sparkle to crystals, etc.  (Only if DISPERSE is defined!)
  32.  *
  33.  *  Arguments for MAT_DIELECTRIC are:
  34.  *    red    grn    blu    rndx    Hartmann
  35.  *
  36.  *  Arguments for MAT_INTERFACE are:
  37.  *    red1    grn1    blu1    rndx1    red2    grn2    blu2    rndx2
  38.  *
  39.  *  The primaries are material transmission per unit length.
  40.  *  MAT_INTERFACE uses dielectric1 for inside and dielectric2 for
  41.  *  outside.
  42.  */
  43.  
  44.  
  45. #define  MLAMBDA    500        /* mean lambda */
  46. #define  MAXLAMBDA    779        /* maximum lambda */
  47. #define  MINLAMBDA    380        /* minimum lambda */
  48.  
  49. #define  MINCOS        0.997        /* minimum dot product for dispersion */
  50.  
  51.  
  52. m_dielectric(m, r)    /* color a ray which hit something transparent */
  53. OBJREC  *m;
  54. register RAY  *r;
  55. {
  56.     double  cos1, cos2, nratio;
  57.     COLOR  mcolor;
  58.     double  mabsorp;
  59.     double  refl, trans;
  60.     FVECT  dnorm;
  61.     double  d1, d2;
  62.     RAY  p;
  63.     register int  i;
  64.  
  65.     if (m->oargs.nfargs != (m->otype==MAT_DIELECTRIC ? 5 : 8))
  66.         objerror(m, USER, "bad arguments");
  67.  
  68.     r->rt = r->rot;                /* just use ray length */
  69.  
  70.     raytexture(r, m->omod);            /* get modifiers */
  71.  
  72.     cos1 = raynormal(dnorm, r);        /* cosine of theta1 */
  73.                         /* index of refraction */
  74.     if (m->otype == MAT_DIELECTRIC)
  75.         nratio = m->oargs.farg[3] + m->oargs.farg[4]/MLAMBDA;
  76.     else
  77.         nratio = m->oargs.farg[3] / m->oargs.farg[7];
  78.     
  79.     if (cos1 < 0.0) {            /* inside */
  80.         cos1 = -cos1;
  81.         dnorm[0] = -dnorm[0];
  82.         dnorm[1] = -dnorm[1];
  83.         dnorm[2] = -dnorm[2];
  84.         setcolor(mcolor, pow(m->oargs.farg[0], r->rot),
  85.                  pow(m->oargs.farg[1], r->rot),
  86.                  pow(m->oargs.farg[2], r->rot));
  87.     } else {                /* outside */
  88.         nratio = 1.0 / nratio;
  89.         if (m->otype == MAT_INTERFACE)
  90.             setcolor(mcolor, pow(m->oargs.farg[4], r->rot),
  91.                      pow(m->oargs.farg[5], r->rot),
  92.                      pow(m->oargs.farg[6], r->rot));
  93.         else
  94.             setcolor(mcolor, 1.0, 1.0, 1.0);
  95.     }
  96.     mabsorp = bright(mcolor);
  97.  
  98.     d2 = 1.0 - nratio*nratio*(1.0 - cos1*cos1);    /* compute cos theta2 */
  99.  
  100.     if (d2 < FTINY)            /* total reflection */
  101.  
  102.         refl = 1.0;
  103.  
  104.     else {                /* refraction occurs */
  105.                     /* compute Fresnel's equations */
  106.         cos2 = sqrt(d2);
  107.         d1 = cos1;
  108.         d2 = nratio*cos2;
  109.         d1 = (d1 - d2) / (d1 + d2);
  110.         refl = d1 * d1;
  111.  
  112.         d1 = 1.0 / cos1;
  113.         d2 = nratio / cos2;
  114.         d1 = (d1 - d2) / (d1 + d2);
  115.         refl += d1 * d1;
  116.  
  117.         refl /= 2.0;
  118.         trans = 1.0 - refl;
  119.  
  120.         if (rayorigin(&p, r, REFRACTED, mabsorp*trans) == 0) {
  121.  
  122.                         /* compute refracted ray */
  123.             d1 = nratio*cos1 - cos2;
  124.             for (i = 0; i < 3; i++)
  125.                 p.rdir[i] = nratio*r->rdir[i] + d1*dnorm[i];
  126.  
  127. #ifdef  DISPERSE
  128.             if (m->otype != MAT_DIELECTRIC
  129.                     || r->rod > 0.0
  130.                     || r->crtype & SHADOW
  131.                     || directinvis
  132.                     || m->oargs.farg[4] == 0.0
  133.                     || !disperse(m, r, p.rdir, trans))
  134. #endif
  135.             {
  136.                 rayvalue(&p);
  137.                 multcolor(mcolor, r->pcol);    /* modify */
  138.                 scalecolor(p.rcol, trans);
  139.                 addcolor(r->rcol, p.rcol);
  140.             }
  141.         }
  142.     }
  143.         
  144.     if (!(r->crtype & SHADOW) &&
  145.             rayorigin(&p, r, REFLECTED, mabsorp*refl) == 0) {
  146.  
  147.                     /* compute reflected ray */
  148.         for (i = 0; i < 3; i++)
  149.             p.rdir[i] = r->rdir[i] + 2.0*cos1*dnorm[i];
  150.  
  151.         rayvalue(&p);            /* reflected ray value */
  152.  
  153.         scalecolor(p.rcol, refl);    /* color contribution */
  154.         addcolor(r->rcol, p.rcol);
  155.     }
  156.  
  157.     multcolor(r->rcol, mcolor);        /* multiply by transmittance */
  158. }
  159.  
  160.  
  161. #ifdef  DISPERSE
  162.  
  163. static
  164. disperse(m, r, vt, tr)        /* check light sources for dispersion */
  165. OBJREC  *m;
  166. RAY  *r;
  167. FVECT  vt;
  168. double  tr;
  169. {
  170.     RAY  sray, *entray;
  171.     FVECT  v1, v2, n1, n2;
  172.     FVECT  dv, v2Xdv;
  173.     double  v2Xdvv2Xdv;
  174.     int  success = 0;
  175.     SRCINDEX  si;
  176.     FVECT  vtmp1, vtmp2;
  177.     double  dtmp1, dtmp2;
  178.     int  l1, l2;
  179.     COLOR  ctmp;
  180.     int  i;
  181.     
  182.     /*
  183.      *     This routine computes dispersion to the first order using
  184.      *  the following assumptions:
  185.      *
  186.      *    1) The dependency of the index of refraction on wavelength
  187.      *        is approximated by Hartmann's equation with lambda0
  188.      *        equal to zero.
  189.      *    2) The entry and exit locations are constant with respect
  190.      *        to dispersion.
  191.      *
  192.      *     The second assumption permits us to model dispersion without
  193.      *  having to sample refracted directions.  We assume that the
  194.      *  geometry inside the material is constant, and concern ourselves
  195.      *  only with the relationship between the entering and exiting ray.
  196.      *  We compute the first derivatives of the entering and exiting
  197.      *  refraction with respect to the index of refraction.  This 
  198.      *  is then used in a first order Taylor series to determine the
  199.      *  index of refraction necessary to send the exiting ray to each
  200.      *  light source.
  201.      *     If an exiting ray hits a light source within the refraction
  202.      *  boundaries, we sum all the frequencies over the disc of the
  203.      *  light source to determine the resulting color.  A smaller light
  204.      *  source will therefore exhibit a sharper spectrum.
  205.      */
  206.  
  207.     if (!(r->crtype & REFRACTED)) {        /* ray started in material */
  208.         VCOPY(v1, r->rdir);
  209.         n1[0] = -r->rdir[0]; n1[1] = -r->rdir[1]; n1[2] = -r->rdir[2];
  210.     } else {
  211.                         /* find entry point */
  212.         for (entray = r; entray->rtype != REFRACTED;
  213.                 entray = entray->parent)
  214.             ;
  215.         entray = entray->parent;
  216.         if (entray->crtype & REFRACTED)        /* too difficult */
  217.             return(0);
  218.         VCOPY(v1, entray->rdir);
  219.         VCOPY(n1, entray->ron);
  220.     }
  221.     VCOPY(v2, vt);            /* exiting ray */
  222.     VCOPY(n2, r->ron);
  223.  
  224.                     /* first order dispersion approx. */
  225.     dtmp1 = DOT(n1, v1);
  226.     dtmp2 = DOT(n2, v2);
  227.     for (i = 0; i < 3; i++)
  228.         dv[i] = v1[i] + v2[i] - n1[i]/dtmp1 - n2[i]/dtmp2;
  229.         
  230.     if (DOT(dv, dv) <= FTINY)    /* null effect */
  231.         return(0);
  232.                     /* compute plane normal */
  233.     fcross(v2Xdv, v2, dv);
  234.     v2Xdvv2Xdv = DOT(v2Xdv, v2Xdv);
  235.  
  236.                     /* check sources */
  237.     initsrcindex(&si);
  238.     while (srcray(&sray, r, &si)) {
  239.  
  240.         if (DOT(sray.rdir, v2) < MINCOS)
  241.             continue;            /* bad source */
  242.                         /* adjust source ray */
  243.  
  244.         dtmp1 = DOT(v2Xdv, sray.rdir) / v2Xdvv2Xdv;
  245.         sray.rdir[0] -= dtmp1 * v2Xdv[0];
  246.         sray.rdir[1] -= dtmp1 * v2Xdv[1];
  247.         sray.rdir[2] -= dtmp1 * v2Xdv[2];
  248.  
  249.         l1 = lambda(m, v2, dv, sray.rdir);    /* mean lambda */
  250.  
  251.         if (l1 > MAXLAMBDA || l1 < MINLAMBDA)    /* not visible */
  252.             continue;
  253.                         /* trace source ray */
  254.         normalize(sray.rdir);
  255.         rayvalue(&sray);
  256.         if (bright(sray.rcol) <= FTINY)    /* missed it */
  257.             continue;
  258.         
  259.         /*
  260.          *    Compute spectral sum over diameter of source.
  261.          *  First find directions for rays going to opposite
  262.          *  sides of source, then compute wavelengths for each.
  263.          */
  264.          
  265.         fcross(vtmp1, v2Xdv, sray.rdir);
  266.         dtmp1 = sqrt(si.dom  / v2Xdvv2Xdv / PI);
  267.  
  268.                             /* compute first ray */
  269.         for (i = 0; i < 3; i++)
  270.             vtmp2[i] = sray.rdir[i] + dtmp1*vtmp1[i];
  271.  
  272.         l1 = lambda(m, v2, dv, vtmp2);        /* first lambda */
  273.         if (l1 < 0)
  274.             continue;
  275.                             /* compute second ray */
  276.         for (i = 0; i < 3; i++)
  277.             vtmp2[i] = sray.rdir[i] - dtmp1*vtmp1[i];
  278.  
  279.         l2 = lambda(m, v2, dv, vtmp2);        /* second lambda */
  280.         if (l2 < 0)
  281.             continue;
  282.                     /* compute color from spectrum */
  283.         if (l1 < l2)
  284.             spec_rgb(ctmp, l1, l2);
  285.         else
  286.             spec_rgb(ctmp, l2, l1);
  287.         multcolor(ctmp, sray.rcol);
  288.         scalecolor(ctmp, tr);
  289.         addcolor(r->rcol, ctmp);
  290.         success++;
  291.     }
  292.     return(success);
  293. }
  294.  
  295.  
  296. static int
  297. lambda(m, v2, dv, lr)            /* compute lambda for material */
  298. register OBJREC  *m;
  299. FVECT  v2, dv, lr;
  300. {
  301.     FVECT  lrXdv, v2Xlr;
  302.     double  dtmp, denom;
  303.     int  i;
  304.  
  305.     fcross(lrXdv, lr, dv);
  306.     for (i = 0; i < 3; i++)    
  307.         if (lrXdv[i] > FTINY || lrXdv[i] < -FTINY)
  308.             break;
  309.     if (i >= 3)
  310.         return(-1);
  311.  
  312.     fcross(v2Xlr, v2, lr);
  313.  
  314.     dtmp = m->oargs.farg[4] / MLAMBDA;
  315.     denom = dtmp + v2Xlr[i]/lrXdv[i] * (m->oargs.farg[3] + dtmp);
  316.  
  317.     if (denom < FTINY)
  318.         return(-1);
  319.  
  320.     return(m->oargs.farg[4] / denom);
  321. }
  322.  
  323. #endif  /* DISPERSE */
  324.