#ifndef WARPING_HPP
#define WARPING_HPP
#include <Eigen/Dense>
#include <cmath>
#include <random>
inline Eigen::Vector3f uniform_sphere(const Eigen::Vector2f& sample){
    using std::acos;
    using std::cos;
    using std::sin;
    std::uniform_real_distribution<float> dis(0, 1);
    float theta = 2 * M_PI * sample(0);
    float phi = acos(2 * sample(1) - 1);
    Eigen::Vector3f ret;
    ret(0) = sin(phi) * cos(theta);
    ret(1) = sin(phi) * sin(theta);
    ret(2) = cos(phi);
    return ret;
}
inline Eigen::Vector3f cosine_hemisphere(const Eigen::Vector3f& normal, const Eigen::Vector2f& sample) {
    using std::asin;
    using std::cos;
    using std::sin;
    using std::sqrt;
	//constexpr float r = 1;
	//float theta = asin(sqrt(sample.y()));
	//float phi = 2 * M_PI * sample.x();
	//float z = r * cos(theta);
	//float _r = sqrt(r * r - z * z);
	//float x = _r * cos(phi);
	//float y = _r * sin(phi);
    float r = sqrt(sample.x());
    float phi = 2 * M_PI * sample.y();
    float x = r * cos(phi);
    float y = r * sin(phi);
    float z = sqrt(1 - r * r);
    Eigen::Vector3f notnormal = normal;
    if(normal.y() == 0){
        std::swap(notnormal.x(), notnormal.z());
        notnormal.z() *= -1;
    }
    else{
        std::swap(notnormal.x(), notnormal.y());
        notnormal.y() *= -1;
    }
    Eigen::Vector3f notnormal2 = normal.cross(notnormal);
    notnormal = notnormal2.cross(normal);
    Eigen::Matrix3f trf;
    trf.col(0) = notnormal;
    trf.col(1) = notnormal2;
    trf.col(2) = normal;
    Eigen::Vector3f ret(x, y, z);
	return Eigen::Vector3f(trf * ret);
}
inline float beckmann_eval(const Eigen::Vector3f& m, float alpha) {
    using std::abs;
    using std::sqrt;
    using std::exp;
    using std::pow;

	if (m.z() > 0 && abs(m.squaredNorm() - 1) < 0.001f) {
		float costheta = m.z();
		float r = sqrt(m.x() * m.x() + m.y() * m.y());
		float tantheta = r / m.z();
		float upper = exp((-tantheta * tantheta) / (alpha * alpha));
		float lower = M_PI * alpha * alpha * pow(costheta, 3);
		return upper / lower;
	}
	return 0.0f;
}

inline Eigen::Vector3f beckmann(const Eigen::Vector2f& sample, float alpha) {
	using std::cos;
	using std::sin;
	using std::sqrt;
    using std::atan;
    using std::log;

	float theta = atan(sqrt(-alpha * alpha * log(1 - sample.y())));
	float pha = 2 * M_PI * sample.x();
	float z = cos(theta);
	float r = sqrt(1 - z * z);
	float y = r * cos(pha);
	float x = r * sin(pha);
	return Eigen::Vector3f(x, y, z);
}

#endif