home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 8 / CDASC08.ISO / NEWS / RADIANCE / LIB / HE.CAL < prev    next >
Text File  |  1993-10-07  |  4KB  |  113 lines
  1. {
  2.     He-Torrance Reflectance Model
  3.  
  4.     This is the simplified version that doesn't account for
  5.   changes in reflection due to changes in wavelength.  Also,
  6.   specular and directional-diffuse hightlights are left uncolored
  7.   because coloring them requires multiple evaluations of some
  8.   very expensive functions.
  9.  
  10.  *  Arguments for the BRDTfunc type are:
  11.  *    10+    rrefl    grefl    brefl
  12.  *        rtrns    gtrns    btrns
  13.  *        rbrtd    gbrtd    bbrtd
  14.  *        funcfile    transform
  15.  *    0
  16.  *    6+    red    grn    blu    rspec    trans    tspec    A7 ..
  17.  *
  18.  *    In addition to the normal variables available to functions,
  19.  *  we can use the following:
  20.  *        NxP, NyP, NzP -        perturbed surface normal
  21.  *        RdotP -            perturbed ray dot product
  22.  *        CrP, CgP, CbP -        perturbed material color
  23. }
  24.  
  25.             { Constants }
  26. lambda : .5;            { wavelength (microns) }
  27. z0err : .0001;        { accepted error in value of z0 }
  28. Dsumlim : .000001;        { last term of D summation }
  29. Dsummax : 200;        { maximum terms in D summation }
  30.  
  31.             { Parameters }
  32. sigma0 = arg(7);        { surface height deviation (microns) }
  33. tau = arg(8);        { correlation distance (microns) }
  34. n_real = arg(9);        { real part of index of refraction }
  35. n_imag = arg(10);        { imaginary part of index of refraction }
  36.             { Derived parameters }
  37. n_k = n_imag/n_real;
  38.  
  39.             { Repeated formulas }
  40. cotexp(t) = tau/sigma0/2/tan(t);
  41. shadowf2(et,erfcet) = (1-.5*erfcet) /
  42.         ((Exp(-sq(et))/sqrt(PI)/et - erfcet)/2 + 1);
  43. shadowf1(t) = or(FTINY-sigma0, .01-abs(t));
  44. shadowf0(t) = abs(t) - (PI/2-.0001);
  45. shadowf(t) = if(shadowf0(t), 0, if(shadowf1(t), 1,
  46.         shadowf2(cotexp(t), erfc(cotexp(t)))));
  47. K(t) = if(abs(t)-FTINY, tan(t) * erfc(cotexp(t)), 0);
  48. fuvA(ct) = sq(n_real)*(1-sq(n_k)) - (1-sq(ct));
  49. fuvB(ct) = sqrt(sq(fuvA(ct)) + 4*sq(sq(n_real)*n_k));
  50. fu2(ct) = (fuvA(ct) + fuvB(ct))/2;
  51. fv2(ct) = (-fuvA(ct) + fuvB(ct))/2;
  52. fperp2(ct) = (sq(ct-sqrt(fu2(ct))) + fv2(ct)) /
  53.         (sq(ct+sqrt(fu2(ct))) + fv2(ct));
  54. fpara2(ct) = (sq(sq(n_real)*(1-sq(n_k))*ct - sqrt(fu2(ct))) +
  55.         sq(2*sq(n_real)*n_k*ct - sqrt(fv2(ct)))) /
  56.         (sq(sq(n_real)*(1-sq(n_k))*ct + sqrt(fu2(ct))) +
  57.         sq(2*sq(n_real)*n_k*ct + sqrt(fv2(ct))));
  58. fresnel2(ct) = (fperp2(ct) + fpara2(ct))/2;
  59.  
  60.             { Formulas dependent only on reflected direction }
  61. theta_r = Acos(RdotP);
  62. shadowf_r = shadowf(theta_r);
  63. K_r = K(theta_r);
  64. srx = Dy*NzP - Dz*NyP; sry = Dz*NxP - Dx*NzP; srz = Dx*NyP - Dy*NxP;
  65. srn2 = sq(srx) + sq(sry) + sq(srz);
  66. prx = sry*Dz - srz*Dy;
  67. pry = srz*Dx - srx*Dz;
  68. prz = srx*Dy - sry*Dx;
  69. s = fresnel2(RdotP)*Exp(-g(RdotP))*sq(shadowf_r);
  70.  
  71.             { Formulas dependent on incident direction }
  72.         { z0 }
  73. z0d(Ki,z) = -(Ki+K_r)/(4*sigma0)*z*Exp(-sq(z/sigma0)/2) - sqrt(PI/2);
  74. z0lim(x) = if(x, max(x,z0err), min(x,-z0err));
  75. z0off(Ki,z) = (sigma0/4*(Ki+K_r)*Exp(-sq(z/sigma0)/2)-sqrt(PI/2)*z)/
  76.         z0lim(z0d(Ki,z));
  77. z0root(Ki, x0, x1, i) = if(i,
  78.             if(z0err-abs(x1-x0),
  79.                 x1,
  80.                 z0root(Ki,x1,x1-z0off(Ki,x1),i-1)),
  81.             0);
  82. z0(ti) = z0root(K(ti), .1, -z0off(K(ti),.1), 100);
  83.         { sigma }
  84. sigma(ti) = if( FTINY-sigma0, sigma0,
  85.         sigma0/sqrt(1+sq(z0(ti)/sigma0)) );
  86.         { g }
  87. g(cti) = sq(2*PI/lambda*sigma(Acos(cti))*(cti+RdotP));
  88.         { |F|^2 }
  89. fresnel2dd(kix,kiy,kiz) = fresnel2(sqrt(sq(kix-Dx) + sq(kiy-Dy) +
  90.                     sq(kiz-Dz))/2);
  91.         { G }
  92. G2( kix,kiy,kiz, six,siy,siz ) =
  93.     sq((sq(Dx-kix)+sq(Dy-kiy)+sq(Dz-kiz))/(Dz-kiz)) /
  94.     sq(sq(Dy*kiz-Dz*kiy)+sq(Dz*kix-Dx*kiz)+sq(Dx*kiy-Dy*kix)) *
  95.     (sq(srx*kix+sry*kiy+srz*kiz) + 
  96.      sq(prx*kix+pry*kiy+prz*kiz)) *
  97.     (sq(six*Dx+siy*Dy+siz*Dz) +
  98.      sq((siy*kiz-siz*kiy)*Dx+(siz*kix-six*kiz)*Dy+(six*kiy-siy*kix)*Dz)) /
  99.     srn2 / (sq(six)+sq(siy)+sq(siz));
  100. G(kix,kiy,kiz) = G2(kix,kiy,kiz,
  101.             kiy*NzP-kiz*NyP, kiz*NxP-kix*NzP, kix*NyP-kiy*NxP);
  102.         { D }
  103. Dsum2(m,lt,c,t,e,g) = if(or(m-Dsummax,and(lt-t,Dsumlim-t)),0,
  104.     t+Dsum2(m+1,t,c*g/(m+1),c*g/(m+1)*Exp(-g-e/(m+1))/(m+1),e,g));
  105. Dsum(e,g) = Dsum2(1,0,g,g*Exp(-g-e),e,g);
  106. D(kix,kiy,kiz) = sq(PI)/4/sq(lambda)*sq(tau) *
  107.     Dsum(sq(2*PI/lambda)/4*sq(tau)*(sq(kix-Dx)+sq(kiy-Dy)),
  108.         g(kix*NxP+kiy*NyP+kiz*NzP));
  109.         { rho_dd }
  110. dd2(cti) = shadowf_r*shadowf(Acos(cti))/cti/RdotP;
  111. dd(kix,kiy,kiz) = dd2(kix*NxP+kiy*NyP+kiz*NzP)*G(kix,kiy,kiz)*
  112.         fresnel2dd(kix,kiy,kiz)/PI*D(kix,kiy,kiz);
  113.