/* ****************************************************************
 *                            F.c
 * ****************************************************************
 *  MODULE PURPOSE:
 *      This module contains routines to for fresnel calculations.
 *
 *  MODULE CONTENTS:
 *      F_d_pl          - reflectance for a dielectric,
 *                          parallel polarized,
 *      F_d_pp          - reflectance for a dielectric,
 *                          perpendicular polarized
 *      F_d_R           - reflectance for a dielectric
 *      F_c_pl          - reflectance for a conductor,
 *                          parallel polarized,
 *      F_c_pp          - reflectance for a conductor,
 *                          perpendicular polarized
 *      F_c_R           - reflectance for a conductor
 *      F_approx_n      - approximate the average index of
 *                          refraction for a material
 *      F_approx_k      - approximate the average absoption
 *                          index for a material
 *      F_approx_Fr     - approximate conductor reflectance
 *      F_approx_Fr_Ft  - approximate dielectric reflectance
 *                          and transmittance
 *
 *
 *  ASSUMPTIONS:
 *      The following are defined in math.h
 *          HUGE    largest floating pont value
 *          M_PI_2  pi/2
 *          M_PI_4  pi/4
 */
#include <stdio.h>
#include <math.h>
#include "geo.h"
#include "F.h"

/* ****************************************************************
 *  double F_d_pl (N, L, T, ni, nt)
 *   DIR_VECT *N, *L, *T    (in) - N, L, T vectors
 *   double   ni        (in) - incident index of refraction
 *   double   nt        (in) - transmitted index of refraction
 *
 *  Returns the reflectance of a dielectric for light with
 *   polarization parallel to the plane of incidence (plane of the
 *   N and L vector), refer to Eq.(\*(eJ). N and L are assumed to
 *   be to the same side of the surface.
 */
double F_d_pl (N, L, T, ni, nt)
DIR_VECT    *N, *L, *T;
double      ni, nt;
{   double  amplitude, N_dot_L, N_dot_T;
    N_dot_L = geo_dot(N, L);
    N_dot_T = geo_dot(N, T);
    amplitude = ((nt * N_dot_L) + (ni * N_dot_T)) /
		((nt * N_dot_L) - (ni * N_dot_T));
    return amplitude * amplitude;
}

/* ****************************************************************
 *  double F_d_pp (N, L, T, ni, nt)
 *   DIR_VECT *N, *L, *T    (in) - N, L, T vectors
 *   double   ni        (in) - incident index of refraction
 *   double   nt        (in) - transmitted index of refraction
 *
 *  Returns the reflectance of a dielectric for light with
 *   polarization perpendicular to the plane of incidence
 *   (plane of the N and L vector), refer to Eq.(\*(eK). N and L
 *   are assumed to be to the same side of the surface.
 */
double F_d_pp (N, L, T, ni, nt)
DIR_VECT    *N, *L, *T;
double      ni, nt;
{   double  amplitude, N_dot_L, N_dot_T;
    N_dot_L = geo_dot(N, L);
    N_dot_T = geo_dot(N, T);
    amplitude = ((ni * N_dot_L) + (nt * N_dot_T)) /
		((ni * N_dot_L) - (nt * N_dot_T));
    return amplitude * amplitude;
}

/* ****************************************************************
 *  double F_d_R (N, L, T, ni, nt)
 *   DIR_VECT *N, *L, *T    (in) - N, L, T vectors
 *   double   ni        (in) incident index of refraction
 *   double   nt        (in) transmitted index of refraction
 *
 *  Returns the average reflectance for a dielectric.  Refer to
 *   Eq.(\*(eL).  N and L are assumed to be to the same side of the
 *   surface.
 */
double F_d_R (N, L, T, ni, nt)
DIR_VECT    *N, *L, *T;
double      ni, nt;
{   return (F_d_pl (N, L, T, ni, nt) +
	F_d_pp (N, L, T, ni, nt)) / 2.0;
}

/* ****************************************************************
 *  double F_c_pl (N, L, nt, k)
 *   DIR_VECT   *N, *L  (in) - N, L vectors
 *   double     nt      (in) - index of refraction
 *   double     k       (in) - absorption coefficient
 *
 *  Returns the reflectance of a conductor for light with
 *   polarization parallel to the plane of incidence (plane of the
 *   N and L vector), refer to Eq.(\*(eO).  N and L are assumed to
 *   be to the same side of the surface.
 */
double F_c_pl (N, L, nt, k)
DIR_VECT    *N, *L;
double      nt, k;
{   double  tmp_f, N_dot_L;
    N_dot_L = geo_dot(N, L);
    tmp_f = ((nt * nt) + (k * k)) * N_dot_L * N_dot_L;
    return (tmp_f - (2.0 * nt * N_dot_L) + 1) /
	   (tmp_f + (2.0 * nt * N_dot_L) + 1);
}

/* ****************************************************************
 *  double F_c_pp (N, L, nt, k)
 *   DIR_VECT   *N, *L  (in) - N, L vectors
 *   double     nt      (in) - index of refraction
 *   double     k       (in) - absorption coefficient
 *
 *  Returns the reflectance of a conductor for light with
 *   polarization perpendicular to the plane of incidence (plane of
 *   the N and L vector), refer to Eq.(\*(eP). N and L are assumed
 *   to be to the same side of the surface.
 */
double F_c_pp (N, L, nt, k)
DIR_VECT    *N, *L;
double      nt, k;
{   double  tmp_f, N_dot_L;
    N_dot_L = geo_dot(N, L);
    tmp_f = (nt * nt) + (k * k);
    return (tmp_f - (2.0 * nt * N_dot_L) + (N_dot_L * N_dot_L)) /
	   (tmp_f + (2.0 * nt * N_dot_L) + (N_dot_L * N_dot_L));
}

/* ****************************************************************
 *  double F_c_R (N, L, nt, k)
 *   DIR_VECT   *N, *L  (in) - N, L vectors
 *   double     nt      (in) - index of refraction
 *   double     k       (in) - absorption coefficient
 *
 *  Returns the average reflectance for a conductor.  Refer to
 *   Eq.(\*(eL).  N and L are assumed to be to the same side of
 *   the surface.
 */
double F_c_R (N, L, nt, k)
DIR_VECT    *N, *L;
double      nt, k;
{   return (F_c_pl (N, L, nt, k) + F_c_pp (N, L, nt, k))/2.0;
}

/* ****************************************************************
 *  double F_approx_n (Ro, n_samp, mtl)
 *   double     *Ro     (out) - average reflectance
 *   int        n_samp  (in)  - number of material samples
 *   double     *mtl    (in)  - material spectral curve
 *
 *  Returns the approximate n, and loads Ro.  This gives the
 *   correct n for a single mtl reflectance if the material is a
 *   dielectric, and the correct n assuming k=0 for a conductor,
 *   refer to Eq.(\*(eN).
 */
double F_approx_n (Ro, n_samp, mtl)
double      *Ro;
int         n_samp;
double      *mtl;
{   int     ct;
    *Ro = 0.0;
    for (ct=n_samp; --ct>=0; ) *Ro += *mtl++;
    *Ro /= (double)n_samp;
    return (1.0 + sqrt(*Ro)) / (1.0 - sqrt(*Ro));
}

/* ****************************************************************
 *  double F_approx_k (Ro, n_samp, mtl)
 *   double     *Ro     (out) - average reflectance
 *   int        n_samp  (in)  - number of material samples
 *   double     *mtl    (in)  - material spectral curve
 *
 *  Returns the approximate k, and loads Ro.  This assumes n=1.0
 *   and is applicable to conductors only. Refer to Eq.(\*(eQ).
 */
double F_approx_k (Ro, n_samp, mtl)
double      *Ro;
int         n_samp;
double      *mtl;
{   int     ct;
    *Ro = 0.0;
    for (ct=n_samp; --ct>=0; ) *Ro += *mtl++;
    *Ro /= (double)n_samp;
    return 2.0 * sqrt(*Ro / (1.0 - *Ro));
}

/* ****************************************************************
 *  double *F_approx_Fr (N, L, Ro, n, k, n_samp, mtl, Fr)
 *   DIR_VECT   *N, *L  (in) - N, L vectors
 *   double     Ro      (in) - average reflectance
 *   double     n       (in) - measured or approximate index
 *                              of refraction
 *   double     k       (in) - measured or approximate
 *                              absorption coefficient
 *   int        n_samp  (in)  number of material samples
 *   double     *mtl    (in)  material spectral curve
 *   double     *Fr     (out) reflectance
 *
 *  Loads the vector Fr with the approximate reflectance for each
 *   color sample point for a conductor, refer to Eq.(\*(eR).
 *   Returns the pointer to Fr.  The measured values for n and k
 *   should be used if they are available.  Otherwise, assume k=0
 *   and approximate n using F_approx_n(), or assume n=1 and
 *   approximate k using F_approx_k().
 */
double *F_approx_Fr (N, L, Ro, n, k, n_samp, mtl, Fr)
DIR_VECT    *N, *L;
double      Ro, n, k;
int         n_samp;
double      *mtl, *Fr;
{   double      R_theta;    /* ave R for N and L*/
    double      factor, *R_out;
    R_theta = F_c_R (N, L, n, k);
    factor = (R_theta - Ro) / (1.0 - Ro);
    for (R_out=Fr ;--n_samp>=0; Fr++, mtl++)
	if ((*Fr = *mtl + ((1.0 - *mtl) * factor)) < 0.0)
	    *Fr = 0.0;
    return R_out;
}

/* ****************************************************************
 *  double *F_approx_Fr_Ft (N, L, T, Ro, ni, nt,
 *                n_samp, *mtl_p, Fr, Ft)
 *  DIR_VECT    *N, *L, *T  (in) - N, L vectors
 *   double     Ro      (in)  average reflectance
 *   double     ni      (in)  ave n for incident material
 *   double     nt      (in)  ave n for transmitted material
 *   int        n_samp  (in)  number of material samples
 *   double     *mtl_p  (in)  pseudo material curve
 *   double     *Fr     (out) reflectance
 *   double     *Tr     (out) transmittance
 *
 *  Loads the vector Fr with the approximate reflectance and Ft
 *   with the approximate transmittance for each sample point for
 *   a conductor, refer to Eq.(\*(eR).  Returns the pointer to Fr.
 *   The measured value for n should be used if it is available,
 *   otherwise use F_approx_n() to approximate the index of
 *   refraction.
 */
double *F_approx_Fr_Ft (N, L, T, Ro, ni, nt, n_samp, mtl_p, Fr, Ft)
DIR_VECT    *N, *L, *T;
double      Ro, ni, nt;
int         n_samp;
double      *mtl_p, *Fr, *Ft;
{   double      R_theta;    /* ave R for N and L*/
    double      factor, *R_out;

    if (T == NULL)      /* internal reflection */
	for (R_out=Fr; --n_samp>=0; *Fr++ = 1.0, *Ft++ = 0.0) ;
    else {
	R_theta = F_d_R (N, L, T, ni, nt);
	factor = (R_theta - Ro) / (1.0 - Ro);
	for (R_out=Fr; --n_samp>=0; mtl_p++) {
	    if ((*Fr = *mtl_p + ((1.0 - *mtl_p) * factor)) < 0.0)
		*Fr = 0.0;
	    *Ft++ = 1.0 - *Fr++;
	}
    }
    return R_out;
}
/* ************************************************************* */