/* ****************************************************************
 *                            G.c
 * ****************************************************************
 *  MODULE PURPOSE:
 *      This module contains geometric attenuation
 *      functions
 *
 *  MODULE CONTENTS:
 *      D_torrance      - function by Torrance and Sparrow
 *      D_Sancer        - function by Sancer
 */
#include <math.h>
#include "geo.h"

#define ROOT_PI     1.7724538509055159

/* ****************************************************************
 * double G_torrance (N, L, V)
 *  DIR_VECT    *N, *L, *V  (in) - N, L, V vectors
 *  double      dummy       (in) - just to make the
 *                              calls identical
 *
 *  Returns geometric attenuation (Torrance and Sparrow
 *   1967), refer to Eq.(\*(eS).  The direction vectors
 *   N, V, and L are assumed to be oriented to the same
 *   side of the surface.
 */
double G_torrance (N, L, V, dummy)
DIR_VECT        *N, *L, *V;
double          dummy;
{
    double      N_dot_H, V_dot_H;
    double      N_dot_V, N_dot_L;
    double      g1, g2;
    DIR_VECT    *H;

    H = geo_H(V,L);
    N_dot_H = geo_dot (N,H);
    V_dot_H = geo_dot (V,H);
    N_dot_V = geo_dot (N,V);
    N_dot_L = geo_dot (N,L);

    if ((g1 = (2.0 * N_dot_H * N_dot_V) /
	V_dot_H) > 1.0) g1 = 1.0;
    if ((g2 = (2.0 * N_dot_H * N_dot_L) /
	V_dot_H) < g1) return g2;
    else return g1;
}



/* ****************************************************************
 * double G_sancer (N, L, V, m_2)
 *  DIR_VECT    *N, *L, *V  (in) - N, L, V vectors
 *  double      m_2         (in) - m squared
 *
 *  Returns geometric attenuation (Sancer 1969), refer to Eq.(\*(rT).
 *   The direction vectors N, V, and L are assumed to be oriented
 *   to the same side of the surface.
 *
 *  NOTE:
 *   > m can be related to beta using D_cook_init()
 *   > This implementation assumes slope relates roughness and
 *      correlation distance as described by Beckmann. This
 *      explains the missing factor of 2 with m_2 when compared
 *      to Eq.(\*(rT).
 */
double G_sancer (N, L, V, m_2)
DIR_VECT        *N, *L, *V;
double          m_2;
{
    double      N_dot_L, N_dot_V;
    double      c1, c2, root_c1, root_c2, ci, cr;

    if (m_2 <= 0.0) return 1.0;
    if ((N_dot_L = geo_dot(N,L)) <= 0.0) return 0.0;
    c1 = (N_dot_L * N_dot_L) / (m_2 * (1.0 - (N_dot_L * N_dot_L)));
    root_c1 = sqrt(c1);
    ci = (exp(-c1) /
	(2.0 * ROOT_PI * root_c1)) - (erfc(root_c1) / 2.0);

    if ((N_dot_V = geo_dot(N,V)) <= 0.0) return 0.0;
    c2 = (N_dot_V * N_dot_V) / (m_2 * (1.0 - (N_dot_V * N_dot_V)));
    root_c2 = sqrt(c2);
    cr = (exp(-c2) /
	(2.0 * ROOT_PI * root_c2)) - (erfc(root_c2) / 2.0);

    return 1.0/(1.0 + ci + cr);
}
/* ************************************************************* */