#pragma once
#include "../../stdafx.h"
#include "../glLibLocals.h"

#include "../glLibError.h"

#include "glLibMathVec.h"

namespace glLib { namespace Math { namespace Random {

template <typename type_float> inline type_float get_uniform(void) {
	return (type_float)(rand())/(type_float)(RAND_MAX);
}
template <typename type_float> inline type_float get_uniform(type_float left, type_float right) {
	return get_uniform<type_float>()*(right-left) + left;
}

/*inline double get_gaussian(void) {
	//http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
	return pow(-2.0*log(get_uniform()),0.5) * cos(2.0*pi*get_uniform());
}
inline double get_gaussian(double mean) {
	return get_gaussian() + mean;
}*/
template <typename type_float> type_float get_gaussian(void) {
	//http://www.taygeta.com/random/gaussian.html
	float x1, x2, w;//, y1, y2;

	do {
		x1 = (type_float)(2.0) * get_uniform<type_float>() - (type_float)(1.0);
		x2 = (type_float)(2.0) * get_uniform<type_float>() - (type_float)(1.0);
		w = x1 * x1 + x2 * x2;
	} while ( w >= (type_float)(1.0) );

	w = sqrt( (type_float)(-2.0) * log(w) / w );
	//y1 = x1 * w;
	//y2 = x2 * w;
	return x1 * w;
}
template <typename type_float> type_float get_gaussian(type_float mean, type_float stddev) {
	return stddev*get_gaussian<type_float>() + mean;
}

template <typename type_float> void get_disk(type_float point[2], type_float radius=1) {
	//http://mathworld.wolfram.com/DiskPointPicking.html
	type_float rand_angle = get_uniform<type_float>((type_float)(0.0),(type_float)(2.0)*::glLib::Math::Misc::pi);
	type_float rand_radius_sqrt = sqrt(get_uniform<type_float>((type_float)(0.0),radius*radius));
	point[0] = rand_radius_sqrt * cos(rand_angle);
	point[1] = rand_radius_sqrt * sin(rand_angle);
}

template <typename type_float> inline void get_sphere(type_float point[3]) {
	point[0]=get_gaussian<type_float>(); point[1]=get_gaussian<type_float>(); point[2]=get_gaussian<type_float>();
	::glLib::Math::Vec::normalize3(point);
}
template <typename type_float> inline void get_hemisphere_x(type_float point[3]) {
	get_sphere(point); point[0] = fabs(point[0]);
}
template <typename type_float> inline void get_hemisphere_y(type_float point[3]) {
	get_sphere(point); point[1] = fabs(point[1]);
}
template <typename type_float> inline void get_hemisphere_z(type_float point[3]) {
	get_sphere(point); point[2] = fabs(point[2]);
}
template <typename type_float> inline void get_hemisphere(type_float point[3], const type_float normal[3]) {
	get_sphere(point);
	if (::glLib::Math::Vec::dot3(normal,point)<(type_float)(0.0)) ::glLib::Math::Vec::reflect_point_plane(point,normal);
}

template <typename type_float> void get_hemisphere_cosine(type_float point[3], type_float radius=1) {
	//Global Illumination Compendium, pg. 19
	type_float rand_angle = get_uniform((type_float)(0.0),(type_float)(2.0)*::glLib::Math::Misc::pi);
	type_float rand_radius = get_uniform((type_float)(0.0),radius*radius);
	type_float rand_radius_sqrt = sqrt(rand_radius);
	point[0] = rand_radius_sqrt * cos(rand_angle);
	point[1] = sqrt((type_float)(1.0)-rand_radius); //i.e., sqrt(1.0-rand_radius_sqrt*rand_radius_sqrt)
	point[2] = rand_radius_sqrt * sin(rand_angle);
}
template <typename type_float> void get_hemisphere_cosine(type_float point[3], const type_float normal[3], type_float radius=1) {
	/*type_float azimuth = get_uniform((type_float)(0.0),(type_float)(2.0)*::glLib::Math::Misc::pi);
	type_float zenith = acos(sqrt(get_uniform<type_float>()));
	static type_float[3] up = {(type_float)(0.0),(type_float)(1.0),(type_float)(0.0)};
	type_float angle_to_normal = acos(dot3(normal,up)); //assumes normal is normalized*/
	//get_hemisphere_cosine(point,radius);

	type_float first_axis[3];
	type_float abs_normal[3] = {fabs(normal[0]),fabs(normal[1]),fabs(normal[2])};
	if      (abs_normal[0]<abs_normal[1] && abs_normal[0]<abs_normal[2]) { first_axis[0]=(type_float)(1.0); first_axis[1]=(type_float)(0.0); first_axis[2]=(type_float)(0.0); }
	else if (abs_normal[1]<abs_normal[2]                               ) { first_axis[0]=(type_float)(0.0); first_axis[1]=(type_float)(1.0); first_axis[2]=(type_float)(0.0); }
	else                                                                 { first_axis[0]=(type_float)(0.0); first_axis[1]=(type_float)(0.0); first_axis[2]=(type_float)(1.0); }

	type_float basis_x[3], basis_z[3];
	::glLib::Math::Vec::cross(normal,first_axis,basis_x); ::glLib::Math::Vec::normalize3(basis_x);
	::glLib::Math::Vec::cross(normal,basis_x,basis_z);// ::glLib::Math::Vec::normalize3(basis_z);

	type_float disk[2];
	get_disk(disk,radius);
	type_float y = sqrt(radius*radius-disk[0]*disk[0]-disk[1]*disk[1]);

	for (char i=0;i<3;++i) {
		point[i] = disk[0]*basis_x[i] + disk[1]*basis_z[i] + y*normal[i];
	}
}
template <typename type_float, typename type_int> void get_hemisphere_cosine_n(type_float point[3], type_int n, type_float radius=1) {
	//Global Illumination Compendium, pg. 19
	type_float rand_angle = get_uniform((type_float)(0.0),(type_float)(2.0)*(type_float)(::glLib::Math::Misc::pi));
	type_float rand_radius = get_uniform((type_float)(0.0),(type_float)(1.0));

	type_float radical = sqrt(   (type_float)(1.0)   -   pow( rand_radius, (type_float)(2.0)/(n+(type_float)(1.0)) )   );

	point[0] = radius * cos(rand_angle) * radical;
	point[1] = radius * pow(rand_radius,(type_float)(1.0)/(n+(type_float)(1.0)));
	point[2] = radius * sin(rand_angle) * radical;
}

/*class glLibMathRandomDistrib {
	private:
		char _type;
		double _param1, _param2;
	public:
		glLibMathRandomDistrib(const glLibFlagRandom type, double param1=0.0f, double param2=1.0f) {
			if (type==GLLIB_RANDOM_UNIFORM) _type = 1;
			else if (type==GLLIB_RANDOM_GAUSSIAN) _type = 2;

			#ifndef GLLIB_ERROR_MODE_REPORT_ERRORS
			else _type = 3;
			#else
			else if (type==GLLIB_RANDOM_NORMAL) _type = 3;
			else throw GLLIB_ERROR_TYPE_ERROR("Random type not recognized!");
			#endif

			_param1 = param1;
			_param2 = param2;
		}
		double sample(void) {
			switch (_type) {
				case 1: return glLibMathRandom::get_uniform(_param1,_param2); break;
				case 2: return glLibMathRandom::get_gaussian(_param1,_param2); break;
				case 3: return glLibMathRandom::get_gaussian(0.0,1.0); break;
			}
		}
};*/

};};};