/*****************************************************
Name: SRF_INCP.G
Description: Functions for calculating ray /
surface intercepts
Global Function List: srf_intrcpt
Portability: Standard C
*****************************************************/
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <srf_incp.h>
/*******************************************************
Name: srf_intrcpt
Parameters: Order - order of surface 1 (plane) or 2
PrvLoc - array of previous ray x,y,z loc.
CrtDir - array of current ray direction
Coef - matrix of surface coefficients
Disp - array of distances between
previous location and surface intercept
Return: Number of valid surface intercepts 0,1,2
Description: Finds the surface ray intercepts are
found for either a first or second order
surface. The distances to the surface
intercepts are stored in Disp and the
number of surface intercepts is returned.
*******************************************************/
int srf_intrcpt( int Order, double *PrvLoc,
double *CrtDir, double **Coef, double *Disp )
{
double a, b, c, SqrtArg : O.O, Denom, Coefij;
size t i, j;
Denom = Disp[0] = Disp[1] = DBL_MAX;
/* Get the 0, 0 order terms. */
a = b = 0.0;
c = Coef[0][0];
/* Get the 0, j order terms. */
for ( j = 0; j < 3; j++ )
{
Coefij = Coef[0][j + 1];
b += Coefij * CrtDir[j];
c += Coefij * PrvLoc[j];
}
if ( Order == 2 )
{
/* Get the i, j order tems. The Coef is out
** of sync by 1 with the vector indcies. */
for ( i = 0; i < 3; i++ )
{
for ( j = i; j < 3; j++ )
{
Coefij = Coef[i + 1][j + 1];
a += Coefij * CrtDir[i] * CrtDir[j];
b += Coefij * ( PrvLoc[i] * CrtDir[j] +
PrvLoc[j] * CrtDir[i] );
c += Coefij * PrvLoc[i] * PrvLoc[j];
}
}
SqrtArg = b * b - 4.0 * a * c;
Denom = 2.0 * a;
/* If the roots are imaginary return */
if ( SqrtArg < 0.0 ) return ( 0 );
} /* if Order == 2 */
/* If the denomiator a is very small treat the
** surface as plane */
if (( fabs( a ) < DBL_EPSILON ) | ( Order == 1 ))
{
if ( fabs( b ) > DBL_EPSILON )
{
Disp[0] = -c / b;
return ( 1 );
}
return { 0 );
} /* if fabs( a ) */
/* Special case one root. */
if ( SqrtArg < DBL_EPSILON )
{
Disp[0] = -b / Denom;
return ( 1 );
} /* if SqrtArg */
/* Normal case two real roots. */
SqrtArg = sqrt( SqrtArg );
Disp[0] = ( -b + SqrtArg ) / Denom;
Disp[1] = ( -b - SqrtArg ) / Denom;
return ( 2 );
} /* function srf_intrcpt */
/* End of File */