#ifndef STIM_CIRCLE_H
#define STIM_CIRCLE_H

#include <stim/cuda/cudatools/callable.h>
#include <stim/math/plane.h>
#include <stim/math/vector.h>
#include <stim/math/triangle.h>
#include <stim/math/constants.h>
#include <assert.h>
#include <algorithm>
#include <iostream>

namespace stim{

template <typename T>
class circle : plane<T>
{

private:
	
	//stim::vec3<T> Y;
	T R;								//radius of the circle

	/*CUDA_CALLABLE void
	init()
	{
		Y = U.cross(N).norm();
	}*/

public:
	using stim::plane<T>::n;
	using stim::plane<T>::P;
	using stim::plane<T>::N;
	using stim::plane<T>::U;
	using stim::plane<T>::rotate;
	using stim::plane<T>::setU;

	using stim::plane<T>::init;

	///base constructor
	///@param th value of the angle of the starting point from 0 to 360.
	CUDA_CALLABLE
	circle() : plane<T>()
	{
		init();
	}

	///create a circle given a size and position in Z space.
	///@param size: size of the rectangle in ND space.
	///@param z_pos z coordinate of the rectangle.
	CUDA_CALLABLE
	circle(T radius, T z_pos = (T)0) : plane<T>(z_pos)
	{
		//center(stim::vec3<T>(0, 0, z_pos));
		//scale(size);
		//init();
		R = radius;
	}

	///create a rectangle from a center point, normal
	///@param c: x,y,z location of the center.
	///@param n: x,y,z direction of the normal.	
	CUDA_CALLABLE
	circle(vec3<T> c, vec3<T> n = vec3<T>(0,0,1)) : plane<T>(n, c)
	{
		//center(c);
		//normal(n);
		//init();
		R = (T)1;
	}

	/*///create a rectangle from a center point, normal, and size
	///@param c: x,y,z location of the center.
	///@param s: size of the rectangle.
	///@param n: x,y,z direction of the normal.
	CUDA_CALLABLE 
	circle(vec3<T> c, T s, vec3<T> n = vec3<T>(0,0,1)) : plane<T>()
	{
		init();
		center(c);
		rotate(n, U, Y);
		scale(s);
	}
	*/

	///create a rectangle from a center point, normal, and size
	///@param c: x,y,z location of the center.
	///@param s: size of the rectangle.
	///@param n: x,y,z direction of the normal.
	///@param u: x,y,z direction for the zero vector (from where the rotation starts)
	CUDA_CALLABLE
	circle(vec3<T> c, T r, vec3<T> n = vec3<T>(0,0,1), vec3<T> u = vec3<T>(1, 0, 0)) : plane<T>()
	{
		P = c;
		N = n;
		setU(u);
		R = r;
		//init();
		//setU(u);
//		U = u;
		//center(c);
		//normal(n);
		//scale(s);
	}

	///scales the circle by a certain factor
	///@param factor: the factor by which the dimensions of the shape are scaled.
	CUDA_CALLABLE
	void scale(T factor)
	{
		//U *= factor;
		//Y *= factor;
		R *= factor;
	}

	///set the radius of circle to a certain value 
	///@param value: the new radius of the circle
	CUDA_CALLABLE
	void set_R(T value) 
	{
		R = value;
	}

	///sets the normal for the cirlce
	///@param n: x,y,z direction of the normal.
	CUDA_CALLABLE void
	normal(vec3<T> n){
		rotate(n);
	}

	///sets the center of the circle.
	///@param n: x,y,z location of the center.
	CUDA_CALLABLE void
	center(vec3<T> p){
		P = p;
	}

	///boolean comparison
	bool
	operator==(const circle<T> & rhs)
	{
		if(P == rhs.P && U == rhs.U)
			return true;
		else
			return false;
	}

	//returns the point in world space corresponding to the polar coordinate (r, theta)
	CUDA_CALLABLE stim::vec3<T>
		p(T r, T theta) {
		T u = r * cos(theta);				//calculate the coordinates in the planar space defined by the circle
		T v = r * sin(theta);

		vec3<T> V = U.cross(N);				//calculate the orthogonal vector V
		return P + U * u + V * v;			//calculate the cartesian coordinate of the specified point
	}

	//returns the point in world space corresponding to the value theta at radius R
	CUDA_CALLABLE stim::vec3<T>
		p(T theta) {
		return p(R, theta);
	}

	//get the world space value given the polar coordinates (r, theta)

	///get the world space value given the planar coordinates a, b in [0, 1]
	/*CUDA_CALLABLE stim::vec3<T> p(T a, T b)
	{
		stim::vec3<T> result;

		vec3<T> A = this->P - this->U * (T)0.5 - Y * (T)0.5;
		result = A + this->U * a + Y * b;
		return result;
	}*/

	///parenthesis operator returns the world space given rectangular coordinates a and b in [0 1]
	CUDA_CALLABLE stim::vec3<T> operator()(T r, T theta)
	{
		return p(r, theta);
	}

	//parenthesis operator returns the world space coordinate at the edge of the circle given theta
	CUDA_CALLABLE stim::vec3<T> operator()(T theta) {
		return p(theta);
	}

	///returns a vector with the points on the initialized circle.
	///connecting the points results in a circle.
	///@param n: integer for the number of points representing the circle.
	std::vector< stim::vec3<T> > points(unsigned n) {
		std::vector< stim::vec3<T> > result(n);				//initialize a vector of n points

		float dt = stim::TAU / n;
		for (unsigned i = 0; i < n; i++)
			result[i] = p(i * dt);							//calculate a point on the edge of the circle
		return result;
	}	

	///returns a vector with the points on the initialized circle
	///connecting the points results in a circle
	///@param n: integer for the number of points representing the circle
	///the only difference between points and glpoints is that the first point appears twice in the returning lists
	std::vector< stim::vec3<T> > glpoints(unsigned n) {
		std::vector< stim::vec3<T> > result(n + 1);
		float dt = stim::TAU / n;
		for (unsigned i = 0; i < n; i++)
			result[i] = p(i * dt);
		result[n] = p(0);									//close the circle!
		return result;
	}
	
	std::string str() const
	{
		std::stringstream ss;
		ss << "r = "<<R<<"  (P=" << P.str() << ", N=" << N.str() << ", U=" << U.str() << ")";
		return ss.str();
	}

};
}
#endif