source-engine/mathlib/spherical.cpp
FluorescentCIAAfricanAmerican 3bf9df6b27 1
2020-04-22 12:56:21 -04:00

125 lines
3.1 KiB
C++

//========= Copyright Valve Corporation, All rights reserved. ============//
//
// Purpose: spherical math routines
//
//=====================================================================================//
#include <math.h>
#include <float.h> // Needed for FLT_EPSILON
#include "basetypes.h"
#include <memory.h>
#include "tier0/dbg.h"
#include "mathlib/mathlib.h"
#include "mathlib/vector.h"
#include "mathlib/spherical_geometry.h"
// memdbgon must be the last include file in a .cpp file!!!
#include "tier0/memdbgon.h"
float s_flFactorials[]={
1.,
1.,
2.,
6.,
24.,
120.,
720.,
5040.,
40320.,
362880.,
3628800.,
39916800.,
479001600.,
6227020800.,
87178291200.,
1307674368000.,
20922789888000.,
355687428096000.,
6402373705728000.,
121645100408832000.,
2432902008176640000.,
51090942171709440000.,
1124000727777607680000.,
25852016738884976640000.,
620448401733239439360000.,
15511210043330985984000000.,
403291461126605635584000000.,
10888869450418352160768000000.,
304888344611713860501504000000.,
8841761993739701954543616000000.,
265252859812191058636308480000000.,
8222838654177922817725562880000000.,
263130836933693530167218012160000000.,
8683317618811886495518194401280000000.
};
float AssociatedLegendrePolynomial( int nL, int nM, float flX )
{
// evaluate associated legendre polynomial at flX, using recurrence relation
float flPmm = 1.;
if ( nM > 0 )
{
float flSomX2 = sqrt( ( 1 - flX ) * ( 1 + flX ) );
float flFact = 1.;
for( int i = 0 ; i < nM; i++ )
{
flPmm *= -flFact * flSomX2;
flFact += 2.0;
}
}
if ( nL == nM )
return flPmm;
float flPmmp1 = flX * ( 2.0 * nM + 1.0 ) * flPmm;
if ( nL == nM + 1 )
return flPmmp1;
float flPll = 0.;
for( int nLL = nM + 2 ; nLL <= nL; nLL++ )
{
flPll = ( ( 2.0 * nLL - 1.0 ) * flX * flPmmp1 - ( nLL + nM - 1.0 ) * flPmm ) * ( 1.0 / ( nLL - nM ) );
flPmm = flPmmp1;
flPmmp1 = flPll;
}
return flPll;
}
static float SHNormalizationFactor( int nL, int nM )
{
double flTemp = ( ( 2. * nL + 1.0 ) * s_flFactorials[ nL - nM ] )/ ( 4. * M_PI * s_flFactorials[ nL + nM ] );
return sqrt( flTemp );
}
#define SQRT_2 1.414213562373095
FORCEINLINE float SphericalHarmonic( int nL, int nM, float flTheta, float flPhi, float flCosTheta )
{
if ( nM == 0 )
return SHNormalizationFactor( nL, 0 ) * AssociatedLegendrePolynomial( nL, nM, flCosTheta );
if ( nM > 0 )
return SQRT_2 * SHNormalizationFactor( nL, nM ) * cos ( nM * flPhi ) *
AssociatedLegendrePolynomial( nL, nM, flCosTheta );
return
SQRT_2 * SHNormalizationFactor( nL, -nM ) * sin( -nM * flPhi ) * AssociatedLegendrePolynomial( nL, -nM, flCosTheta );
}
float SphericalHarmonic( int nL, int nM, float flTheta, float flPhi )
{
return SphericalHarmonic( nL, nM, flTheta, flPhi, cos( flTheta ) );
}
float SphericalHarmonic( int nL, int nM, Vector const &vecDirection )
{
Assert( fabs( VectorLength( vecDirection ) - 1.0 ) < 0.0001 );
float flPhi = acos( vecDirection.z );
float flTheta = 0;
float S = Square( vecDirection.x ) + Square( vecDirection.y );
if ( S > 0 )
{
flTheta = atan2( vecDirection.y, vecDirection.x );
}
return SphericalHarmonic( nL, nM, flTheta, flPhi, cos( flTheta ) );
}