home *** CD-ROM | disk | FTP | other *** search
/ Compendium Deluxe 2 / LSD and 17bit Compendium Deluxe - Volume II.iso / a / prog / cprog / illum.lha / F.c < prev    next >
Encoding:
C/C++ Source or Header  |  1988-05-18  |  9.9 KB  |  274 lines
  1. /* ****************************************************************
  2.  *                            F.c
  3.  * ****************************************************************
  4.  *  MODULE PURPOSE:
  5.  *      This module contains routines to for fresnel calculations.
  6.  *
  7.  *  MODULE CONTENTS:
  8.  *      F_d_pl          - reflectance for a dielectric,
  9.  *                          parallel polarized,
  10.  *      F_d_pp          - reflectance for a dielectric,
  11.  *                          perpendicular polarized
  12.  *      F_d_R           - reflectance for a dielectric
  13.  *      F_c_pl          - reflectance for a conductor,
  14.  *                          parallel polarized,
  15.  *      F_c_pp          - reflectance for a conductor,
  16.  *                          perpendicular polarized
  17.  *      F_c_R           - reflectance for a conductor
  18.  *      F_approx_n      - approximate the average index of
  19.  *                          refraction for a material
  20.  *      F_approx_k      - approximate the average absoption
  21.  *                          index for a material
  22.  *      F_approx_Fr     - approximate conductor reflectance
  23.  *      F_approx_Fr_Ft  - approximate dielectric reflectance
  24.  *                          and transmittance
  25.  *
  26.  *
  27.  *  ASSUMPTIONS:
  28.  *      The following are defined in math.h
  29.  *          HUGE    largest floating pont value
  30.  *          M_PI_2  pi/2
  31.  *          M_PI_4  pi/4
  32.  */
  33. #include <stdio.h>
  34. #include <math.h>
  35. #include "geo.h"
  36. #include "F.h"
  37.  
  38. /* ****************************************************************
  39.  *  double F_d_pl (N, L, T, ni, nt)
  40.  *   DIR_VECT *N, *L, *T    (in) - N, L, T vectors
  41.  *   double   ni        (in) - incident index of refraction
  42.  *   double   nt        (in) - transmitted index of refraction
  43.  *
  44.  *  Returns the reflectance of a dielectric for light with
  45.  *   polarization parallel to the plane of incidence (plane of the
  46.  *   N and L vector), refer to Eq.(\*(eJ). N and L are assumed to
  47.  *   be to the same side of the surface.
  48.  */
  49. double F_d_pl (N, L, T, ni, nt)
  50. DIR_VECT    *N, *L, *T;
  51. double      ni, nt;
  52. {   double  amplitude, N_dot_L, N_dot_T;
  53.     N_dot_L = geo_dot(N, L);
  54.     N_dot_T = geo_dot(N, T);
  55.     amplitude = ((nt * N_dot_L) + (ni * N_dot_T)) /
  56.         ((nt * N_dot_L) - (ni * N_dot_T));
  57.     return amplitude * amplitude;
  58. }
  59.  
  60. /* ****************************************************************
  61.  *  double F_d_pp (N, L, T, ni, nt)
  62.  *   DIR_VECT *N, *L, *T    (in) - N, L, T vectors
  63.  *   double   ni        (in) - incident index of refraction
  64.  *   double   nt        (in) - transmitted index of refraction
  65.  *
  66.  *  Returns the reflectance of a dielectric for light with
  67.  *   polarization perpendicular to the plane of incidence
  68.  *   (plane of the N and L vector), refer to Eq.(\*(eK). N and L
  69.  *   are assumed to be to the same side of the surface.
  70.  */
  71. double F_d_pp (N, L, T, ni, nt)
  72. DIR_VECT    *N, *L, *T;
  73. double      ni, nt;
  74. {   double  amplitude, N_dot_L, N_dot_T;
  75.     N_dot_L = geo_dot(N, L);
  76.     N_dot_T = geo_dot(N, T);
  77.     amplitude = ((ni * N_dot_L) + (nt * N_dot_T)) /
  78.         ((ni * N_dot_L) - (nt * N_dot_T));
  79.     return amplitude * amplitude;
  80. }
  81.  
  82. /* ****************************************************************
  83.  *  double F_d_R (N, L, T, ni, nt)
  84.  *   DIR_VECT *N, *L, *T    (in) - N, L, T vectors
  85.  *   double   ni        (in) incident index of refraction
  86.  *   double   nt        (in) transmitted index of refraction
  87.  *
  88.  *  Returns the average reflectance for a dielectric.  Refer to
  89.  *   Eq.(\*(eL).  N and L are assumed to be to the same side of the
  90.  *   surface.
  91.  */
  92. double F_d_R (N, L, T, ni, nt)
  93. DIR_VECT    *N, *L, *T;
  94. double      ni, nt;
  95. {   return (F_d_pl (N, L, T, ni, nt) +
  96.     F_d_pp (N, L, T, ni, nt)) / 2.0;
  97. }
  98.  
  99. /* ****************************************************************
  100.  *  double F_c_pl (N, L, nt, k)
  101.  *   DIR_VECT   *N, *L  (in) - N, L vectors
  102.  *   double     nt      (in) - index of refraction
  103.  *   double     k       (in) - absorption coefficient
  104.  *
  105.  *  Returns the reflectance of a conductor for light with
  106.  *   polarization parallel to the plane of incidence (plane of the
  107.  *   N and L vector), refer to Eq.(\*(eO).  N and L are assumed to
  108.  *   be to the same side of the surface.
  109.  */
  110. double F_c_pl (N, L, nt, k)
  111. DIR_VECT    *N, *L;
  112. double      nt, k;
  113. {   double  tmp_f, N_dot_L;
  114.     N_dot_L = geo_dot(N, L);
  115.     tmp_f = ((nt * nt) + (k * k)) * N_dot_L * N_dot_L;
  116.     return (tmp_f - (2.0 * nt * N_dot_L) + 1) /
  117.        (tmp_f + (2.0 * nt * N_dot_L) + 1);
  118. }
  119.  
  120. /* ****************************************************************
  121.  *  double F_c_pp (N, L, nt, k)
  122.  *   DIR_VECT   *N, *L  (in) - N, L vectors
  123.  *   double     nt      (in) - index of refraction
  124.  *   double     k       (in) - absorption coefficient
  125.  *
  126.  *  Returns the reflectance of a conductor for light with
  127.  *   polarization perpendicular to the plane of incidence (plane of
  128.  *   the N and L vector), refer to Eq.(\*(eP). N and L are assumed
  129.  *   to be to the same side of the surface.
  130.  */
  131. double F_c_pp (N, L, nt, k)
  132. DIR_VECT    *N, *L;
  133. double      nt, k;
  134. {   double  tmp_f, N_dot_L;
  135.     N_dot_L = geo_dot(N, L);
  136.     tmp_f = (nt * nt) + (k * k);
  137.     return (tmp_f - (2.0 * nt * N_dot_L) + (N_dot_L * N_dot_L)) /
  138.        (tmp_f + (2.0 * nt * N_dot_L) + (N_dot_L * N_dot_L));
  139. }
  140.  
  141. /* ****************************************************************
  142.  *  double F_c_R (N, L, nt, k)
  143.  *   DIR_VECT   *N, *L  (in) - N, L vectors
  144.  *   double     nt      (in) - index of refraction
  145.  *   double     k       (in) - absorption coefficient
  146.  *
  147.  *  Returns the average reflectance for a conductor.  Refer to
  148.  *   Eq.(\*(eL).  N and L are assumed to be to the same side of
  149.  *   the surface.
  150.  */
  151. double F_c_R (N, L, nt, k)
  152. DIR_VECT    *N, *L;
  153. double      nt, k;
  154. {   return (F_c_pl (N, L, nt, k) + F_c_pp (N, L, nt, k))/2.0;
  155. }
  156.  
  157. /* ****************************************************************
  158.  *  double F_approx_n (Ro, n_samp, mtl)
  159.  *   double     *Ro     (out) - average reflectance
  160.  *   int        n_samp  (in)  - number of material samples
  161.  *   double     *mtl    (in)  - material spectral curve
  162.  *
  163.  *  Returns the approximate n, and loads Ro.  This gives the
  164.  *   correct n for a single mtl reflectance if the material is a
  165.  *   dielectric, and the correct n assuming k=0 for a conductor,
  166.  *   refer to Eq.(\*(eN).
  167.  */
  168. double F_approx_n (Ro, n_samp, mtl)
  169. double      *Ro;
  170. int         n_samp;
  171. double      *mtl;
  172. {   int     ct;
  173.     *Ro = 0.0;
  174.     for (ct=n_samp; --ct>=0; ) *Ro += *mtl++;
  175.     *Ro /= (double)n_samp;
  176.     return (1.0 + sqrt(*Ro)) / (1.0 - sqrt(*Ro));
  177. }
  178.  
  179. /* ****************************************************************
  180.  *  double F_approx_k (Ro, n_samp, mtl)
  181.  *   double     *Ro     (out) - average reflectance
  182.  *   int        n_samp  (in)  - number of material samples
  183.  *   double     *mtl    (in)  - material spectral curve
  184.  *
  185.  *  Returns the approximate k, and loads Ro.  This assumes n=1.0
  186.  *   and is applicable to conductors only. Refer to Eq.(\*(eQ).
  187.  */
  188. double F_approx_k (Ro, n_samp, mtl)
  189. double      *Ro;
  190. int         n_samp;
  191. double      *mtl;
  192. {   int     ct;
  193.     *Ro = 0.0;
  194.     for (ct=n_samp; --ct>=0; ) *Ro += *mtl++;
  195.     *Ro /= (double)n_samp;
  196.     return 2.0 * sqrt(*Ro / (1.0 - *Ro));
  197. }
  198.  
  199. /* ****************************************************************
  200.  *  double *F_approx_Fr (N, L, Ro, n, k, n_samp, mtl, Fr)
  201.  *   DIR_VECT   *N, *L  (in) - N, L vectors
  202.  *   double     Ro      (in) - average reflectance
  203.  *   double     n       (in) - measured or approximate index
  204.  *                              of refraction
  205.  *   double     k       (in) - measured or approximate
  206.  *                              absorption coefficient
  207.  *   int        n_samp  (in)  number of material samples
  208.  *   double     *mtl    (in)  material spectral curve
  209.  *   double     *Fr     (out) reflectance
  210.  *
  211.  *  Loads the vector Fr with the approximate reflectance for each
  212.  *   color sample point for a conductor, refer to Eq.(\*(eR).
  213.  *   Returns the pointer to Fr.  The measured values for n and k
  214.  *   should be used if they are available.  Otherwise, assume k=0
  215.  *   and approximate n using F_approx_n(), or assume n=1 and
  216.  *   approximate k using F_approx_k().
  217.  */
  218. double *F_approx_Fr (N, L, Ro, n, k, n_samp, mtl, Fr)
  219. DIR_VECT    *N, *L;
  220. double      Ro, n, k;
  221. int         n_samp;
  222. double      *mtl, *Fr;
  223. {   double      R_theta;    /* ave R for N and L*/
  224.     double      factor, *R_out;
  225.     R_theta = F_c_R (N, L, n, k);
  226.     factor = (R_theta - Ro) / (1.0 - Ro);
  227.     for (R_out=Fr ;--n_samp>=0; Fr++, mtl++)
  228.     if ((*Fr = *mtl + ((1.0 - *mtl) * factor)) < 0.0)
  229.         *Fr = 0.0;
  230.     return R_out;
  231. }
  232.  
  233. /* ****************************************************************
  234.  *  double *F_approx_Fr_Ft (N, L, T, Ro, ni, nt,
  235.  *                n_samp, *mtl_p, Fr, Ft)
  236.  *  DIR_VECT    *N, *L, *T  (in) - N, L vectors
  237.  *   double     Ro      (in)  average reflectance
  238.  *   double     ni      (in)  ave n for incident material
  239.  *   double     nt      (in)  ave n for transmitted material
  240.  *   int        n_samp  (in)  number of material samples
  241.  *   double     *mtl_p  (in)  pseudo material curve
  242.  *   double     *Fr     (out) reflectance
  243.  *   double     *Tr     (out) transmittance
  244.  *
  245.  *  Loads the vector Fr with the approximate reflectance and Ft
  246.  *   with the approximate transmittance for each sample point for
  247.  *   a conductor, refer to Eq.(\*(eR).  Returns the pointer to Fr.
  248.  *   The measured value for n should be used if it is available,
  249.  *   otherwise use F_approx_n() to approximate the index of
  250.  *   refraction.
  251.  */
  252. double *F_approx_Fr_Ft (N, L, T, Ro, ni, nt, n_samp, mtl_p, Fr, Ft)
  253. DIR_VECT    *N, *L, *T;
  254. double      Ro, ni, nt;
  255. int         n_samp;
  256. double      *mtl_p, *Fr, *Ft;
  257. {   double      R_theta;    /* ave R for N and L*/
  258.     double      factor, *R_out;
  259.  
  260.     if (T == NULL)      /* internal reflection */
  261.     for (R_out=Fr; --n_samp>=0; *Fr++ = 1.0, *Ft++ = 0.0) ;
  262.     else {
  263.     R_theta = F_d_R (N, L, T, ni, nt);
  264.     factor = (R_theta - Ro) / (1.0 - Ro);
  265.     for (R_out=Fr; --n_samp>=0; mtl_p++) {
  266.         if ((*Fr = *mtl_p + ((1.0 - *mtl_p) * factor)) < 0.0)
  267.         *Fr = 0.0;
  268.         *Ft++ = 1.0 - *Fr++;
  269.     }
  270.     }
  271.     return R_out;
  272. }
  273. /* ************************************************************* */
  274.