#ifndef STIM_TENSOR2_H
#define STIM_TENSOR2_H

#include "matrix_sym.h"

namespace stim {

/*This class represents a symmetric rank-2 2D tensor, useful for structure tensors
*/
template<typename T>
class tensor2 : public matrix_sym<T, 2> {

protected:

public:

	//calculate the eigenvectors and eigenvalues of the tensor
	CUDA_CALLABLE void eig(stim::matrix<T, 2>& v, stim::matrix<T, 2>& lambda) {
		
		lambda = 0;							//initialize the eigenvalue matrix to zero
		
		T t = M[0] + M[2];					//calculate the trace of the tensor
		T d = M[0] * M[2] - M[1] * M[1];	//calculate the determinant of the tensor
		
		lambda(0, 0) = t / 2 + sqrt(t*t / 4 - d);
		lambda(1, 1) = t / 2 - sqrt(t*t / 4 - d);

		if (M[1] == 0) {
			v = stim::matrix<T, 2>::identity();
		}
		else {
			v(0, 0) = lambda(0, 0) - d;
			v(0, 1) = lambda(1, 1) - d;
			v(1, 0) = v(1, 1) = M[1];
		}
	}

	CUDA_CALLABLE tensor2<T> operator=(stim::matrix_sym<T, 2> rhs){
		stim::matrix_sym<T, 2>::operator=(rhs);
		return *this;
	}
};


}	//end namespace stim


#endif