#ifndef RTS_PLANE_H
#define RTS_PLANE_H

#include <iostream>
#include <stim/math/vector.h>
#include "rts/cuda/callable.h"


namespace stim{
template <typename T, int D> class plane;
}

template <typename T, int D>
CUDA_CALLABLE stim::plane<T, D> operator-(stim::plane<T, D> v);

namespace stim{

template <class T, int D = 3>
class plane{

	//a plane is defined by a point and a normal

private:

	vec<T, D> P;	//point on the plane
	vec<T, D> N;	//plane normal

	CUDA_CALLABLE void init(){
		P = vec<T, D>(0, 0, 0);
		N = vec<T, D>(0, 0, 1);
	}


public:

	//default constructor
	CUDA_CALLABLE plane(){
		init();
	}
	
	CUDA_CALLABLE plane(vec<T, D> n, vec<T, D> p = vec<T, D>(0, 0, 0)){
		P = p;
		N = n.norm();
	}

	CUDA_CALLABLE plane(T z_pos){
		init();
		P[2] = z_pos;
	}

	//create a plane from three points (a triangle)
	CUDA_CALLABLE plane(vec<T, D> a, vec<T, D> b, vec<T, D> c){
		P = c;
		N = (c - a).cross(b - a);
		if(N.len() == 0)	//handle the degenerate case when two vectors are the same, N = 0
			N = 0;
		else
			N = N.norm();
	}

	template< typename U >
	CUDA_CALLABLE operator plane<U, D>(){

		plane<U, D> result(N, P);
		return result;
	}

	CUDA_CALLABLE vec<T, D> norm(){
		return N;
	}

	CUDA_CALLABLE vec<T, D> p(){
		return P;
	}

	//flip the plane front-to-back
	CUDA_CALLABLE plane<T, D> flip(){
		plane<T, D> result = *this;
		result.N = -result.N;
		return result;
	}

	//determines how a vector v intersects the plane (1 = intersects front, 0 = within plane, -1 = intersects back)
	CUDA_CALLABLE int face(vec<T, D> v){
		
		T dprod = v.dot(N);		//get the dot product between v and N

		//conditional returns the appropriate value
		if(dprod < 0)
			return 1;
		else if(dprod > 0)
			return -1;
		else
			return 0;
	}

	//determine on which side of the plane a point lies (1 = front, 0 = on the plane, -1 = back)
	CUDA_CALLABLE int side(vec<T, D> p){

		vec<T, D> v = p - P;	//get the vector from P to the query point p

		return face(v);
	}

	//compute the component of v that is perpendicular to the plane
	CUDA_CALLABLE vec<T, D> perpendicular(vec<T, D> v){
		return N * v.dot(N);
	}

	//compute the projection of v in the plane
	CUDA_CALLABLE vec<T, D> parallel(vec<T, D> v){
		return v - perpendicular(v);
	}

	CUDA_CALLABLE void decompose(vec<T, D> v, vec<T, D>& para, vec<T, D>& perp){
		perp = N * v.dot(N);
		para = v - perp;
	}

	//get both the parallel and perpendicular components of a vector v w.r.t. the plane
	CUDA_CALLABLE void project(vec<T, D> v, vec<T, D> &v_par, vec<T, D> &v_perp){

		v_perp = v.dot(N);
		v_par = v - v_perp;
	}

	//compute the reflection of v off of the plane
	CUDA_CALLABLE vec<T, D> reflect(vec<T, D> v){

		//compute the reflection using N_prime as the plane normal
		vec<T, D> par = parallel(v);
		vec<T, D> r = (-v) + par * 2;

		/*std::cout<<"----------------REFLECT-----------------------------"<<std::endl;
		std::cout<<str()<<std::endl;
		std::cout<<"v: "<<v<<std::endl;
		std::cout<<"r: "<<r<<std::endl;
		std::cout<<"Perpendicular: "<<perpendicular(v)<<std::endl;
		std::cout<<"Parallel: "<<par<<std::endl;*/
		return r;

	}

	CUDA_CALLABLE rts::plane<T, D> operator-()
	{
		rts::plane<T, D> p = *this;

		//negate the normal vector
		p.N = -p.N;

		return p;
	}

	//output a string
	std::string str(){
		std::stringstream ss;
		ss<<"P: "<<P<<std::endl;
		ss<<"N: "<<N;
		return ss.str();
	}

	///////Friendship
	//friend CUDA_CALLABLE rts::plane<T, D> operator- <> (rts::plane<T, D> v);



};

}

//arithmetic operators

//negative operator flips the plane (front to back)
//template <typename T, int D>




#endif