/*	Copyright (c) 2004 Simon Brown

	Permission is hereby granted, free of charge, to any person obtaining a 
	copy of this software and associated documentation files (the "Software"),
	to deal in the Software without restriction, including without limitation
	the rights to use, copy, modify, merge, publish, distribute, sublicense,
	and/or sell copies of the Software, and to permit persons to whom the
	Software is furnished to do so, subject to the following conditions:

	The above copyright notice and this permission notice shall be included in
	all copies or substantial portions of the Software.

	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
	IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
	FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
	THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 
	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 
	FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 
	DEALINGS IN THE SOFTWARE.
*/

#include <iostream>
#include <cmath>

static double const pi = 3.14159265358979323433832792;

/*inline double harmonic( int n, double phi, double theta )
{
	for( int l = 0;; ++l )
	{
		if( std::abs( n ) <= l )
			return harmonic( l, n, phi, theta );
		n -= ( 2*l + 1 );
	}
}*/
	
inline double harmonic( int l, int m, double phi, double theta )
{
	static double const root_pi = std::sqrt( pi );
	
	switch( l )
	{
	case 0:
		{
			static double const a0 = 1.0 / ( 2.0 * root_pi );
			return a0;
		}
		break;

	case 1:
		{
			static double const a0 = std::sqrt( 3.0 ) / ( 2.0 * root_pi );

			double const ct = std::cos( theta );
			double const st = std::sin( theta );
			double const cp = std::cos( phi );
			double const sp = std::sin( phi );
			
			switch( m )
			{
			case -1:
				return a0*sp*st;
				
			case 0:
				return a0*ct;

			default: // case 1
				return a0*cp*st;
			}
		}
		break;

	case 2:
		{
			static double const a0 = std::sqrt( 15.0 ) / ( 2.0 * root_pi );
			static double const a1 = std::sqrt( 5.0 ) / ( 4.0 * root_pi );
			static double const a2 = std::sqrt( 15.0 ) / ( 4.0 * root_pi );

			double const ct = std::cos( theta );
			double const st = std::sin( theta );
			double const cp = std::cos( phi );
			double const sp = std::sin( phi );
			
			switch( m )
			{
			case -2:
				return a0*cp*sp*st*st;
				
			case -1:
				return a0*sp*ct*st;
				
			case 0:
				return a1*( 3.0*ct*ct - 1.0 );

			case 1:
				return a0*cp*ct*st;
				
			default: // case 2
				return a2*( cp*cp - sp*sp )*st*st;
			}
		}
		break;

	case 3:
		{
			static double const a0 = std::sqrt( 35.0 ) / ( 4.0 * std::sqrt( 2.0 * pi ) );
			static double const a1 = std::sqrt( 105.0 ) / ( 2.0 * root_pi );
			static double const a2 = std::sqrt( 21.0 ) / ( 4.0 * std::sqrt( 2.0 * pi ) );
			static double const a3 = std::sqrt( 7.0 ) / ( 4.0 * root_pi );
			static double const a4 = std::sqrt( 105.0 ) / ( 4.0 * root_pi );
			
			double const ct = std::cos( theta );
			double const st = std::sin( theta );
			double const cp = std::cos( phi );
			double const sp = std::sin( phi );
			
			switch( m )
			{
			case -3:
				return a0*( 3.0*cp*cp - sp*sp )*sp*st*st*st;
				
			case -2:
				return a1*cp*sp*ct*st*st;

			case -1:
				return a2*( 5.0*ct*ct - 1.0 )*sp*st;

			case 0:
				return a3*( 5.0*ct*ct*ct - 3.0*ct );

			case 1:
				return a2*( 5.0*ct*ct - 1.0 )*cp*st;

			case 2:
				return a4*( cp*cp - sp*sp )*ct*st*st;

			default: // case 3
				return a0*( cp*cp - 3.0*sp*sp )*cp*st*st*st;
			}			
		}
		break;

	case 4:
		{
			static double const a0 = ( 3.0 * std::sqrt( 35.0 ) ) / ( 4.0 * root_pi );
			static double const a1 = ( 3.0 * std::sqrt( 35.0 ) ) / ( 4.0 * std::sqrt( 2.0 * pi ) );
			static double const a2 = ( 3.0 * std::sqrt( 5.0 ) ) / ( 4.0 * root_pi );
			static double const a3 = ( 3.0 * std::sqrt( 5.0 ) ) / ( 4.0 * std::sqrt( 2.0 * pi ) );
			static double const a4 = 3.0 / ( 16.0 * root_pi );
			static double const a5 = ( 3.0 * std::sqrt( 5.0 ) ) / ( 8.0 * root_pi );
			static double const a6 = ( 3.0 * std::sqrt( 35.0 ) ) / ( 16.0 * root_pi );

			double const ct = std::cos( theta );
			double const st = std::sin( theta );
			double const cp = std::cos( phi );
			double const sp = std::sin( phi );
			
			switch( m )
			{
			case -4:
				return a0*( cp*cp - sp*sp )*cp*sp*st*st*st*st;

			case -3:
				return a1*( 3.0*cp*cp - sp*sp )*sp*ct*st*st*st;

			case -2:
				return a2*( 7.0*ct*ct - 1.0 )*cp*sp*st*st;

			case -1:
				return a3*( 7.0*ct*ct - 3.0 )*sp*ct*st;

			case 0:
				return a4*( 3.0 - 30.0*ct*ct + 35.0*ct*ct*ct*ct );

			case 1:
				return a3*( 7.0*ct*ct - 3.0 )*cp*ct*st;

			case 2:
				return a5*( 7.0*ct*ct - 1.0 )*( cp*cp - sp*sp )*st*st;

			case 3:
				return a1*( cp*cp - 3.0*sp*sp )*cp*ct*st*st*st;

			default: // case 4
				return a6*( sp*sp*sp*sp - 6.0*cp*cp*sp*sp + cp*cp*cp*cp )*st*st*st*st;
			}
		}
		break;

	default:
		return 0.0;
	}
}

inline double convolve( int l1, int m1, int l2, int m2 )
{
	static const int phi_levels = 72;
	static const int theta_levels = 72;

	double sum( 0.0 );
	for( int theta_count = 0; theta_count < theta_levels; ++theta_count )
	{
		double theta = ( pi * double( theta_count ) ) / double( theta_levels );
		double sin_theta = std::sin( theta );

		for( int phi_count = 0; phi_count < phi_levels; ++phi_count )
		{
			double phi = ( 2.0 * pi * double( phi_count ) ) / double( phi_levels );

			sum += harmonic( l1, m1, phi, theta ) * harmonic( l2, m2, phi, theta ) * sin_theta;
		}
	}

	double areaPerSample = ( 2.0 * pi * pi ) / double( phi_levels * theta_levels );
	return sum * areaPerSample;
}

void test_orthogonality( int max_order )
{
	std::cout << std::endl;
	
	bool passed = true;
	for( int l1 = 0; l1 <= max_order; ++l1 )
	{
		std::cout << "testing order " << l1 << "..." << std::endl;
		for( int m1 = -l1; m1 <= l1; ++m1 )
		{
			for( int l2 = 0; l2 <= max_order; ++l2 )
			{
				for( int m2 = -l2; m2 <= l2; ++m2 )
				{
					double sum = convolve( l1, m1, l2, m2 );
					
					bool zero = std::abs( sum ) < 0.01;
					bool one = std::abs( 1.0 - sum ) < 0.01;
					
					if( ( m1 == m2 && l1 == l2 && !one )
						|| ( ( m1 != m2 || l1 != l2 ) && !zero ) )
					{
						std::cout << "(" << l1 << "," << m1 << ")*"
							<< "(" << l2 << "," << m2 << ") = " << sum << std::endl;
						passed = false;
					}
				}
			}
		}
	}

	std::cout << "orthogonality test: " << ( passed ? "passed" : "failed" ) << std::endl;
}

void test_irradiance( int max_order, double light_phi, double light_theta )
{
	std::cout << std::endl;
	
	static const int phi_levels = 72;
	static const int theta_levels = 72;

	double a[] = { pi, 2.0*pi/3.0, pi/4.0, 0.0, -pi/24.0 };

	double error_sum = 0.0;
	bool passed = true;
	for( int theta_count = 0; theta_count < theta_levels; ++theta_count )
	{
		double sample_theta = ( pi * double( theta_count ) ) / double( theta_levels );
		for( int phi_count = 0; phi_count < phi_levels; ++phi_count )
		{
			double sample_phi = ( 2.0 * pi * double( phi_count ) ) / double( phi_levels );

			double irradiance = 0.0;
			for( int l = 0; l <= max_order; ++l )
			{
				for( int m = -l; m <= l; ++m )
				{
					irradiance += a[l]
						* harmonic( l, m, light_phi, light_theta )
						* harmonic( l, m, sample_phi, sample_theta );
				}
			}

			double lambertian = 
				std::sin( light_theta )*std::cos( light_phi )*std::sin( sample_theta )*std::cos( sample_phi )
				+ std::sin( light_theta )*std::sin( light_phi )*std::sin( sample_theta )*std::sin( sample_phi )
				+ std::cos( light_theta )*std::cos( sample_theta );
			if( lambertian < 0.0 )
				lambertian = 0.0;

			double error = std::abs( irradiance - lambertian );
			error_sum += error;

			if( error > 0.1 )
			{
				std::cout << "irradiance: " << irradiance << ", lambertian: " << lambertian << std::endl;
				passed = false; 
			}
		}
	}
	std::cout << "irradiance average error: " 
		<< ( error_sum / double( theta_levels*phi_levels ) ) << std::endl;

	std::cout << "irradiance test: " << ( passed ? "passed" : "failed" ) << std::endl;
}

int main()
{
	// test orthogonality of the basis
	test_orthogonality( 4 );

	// test a couple of light orientations
	test_irradiance( 2, pi/4.0, pi/4.0 );
	test_irradiance( 2, 0, pi/3.0 );

	return 0;
}