#include"../include/elastic.h"
using std::cout;
using std::endl;
using std::vector;
using std::string;
using namespace elas;
vector<vector<int> > deformation_type{
		{1,2,3,4,5,6},
		{-2,1,4,-3,6,-5},
		{3,-5,-1,6,2,-4},
		{-4,-6,5,1,-3,2},
		{5,4,6,-2,-1,-3},
		{-6,3,-2,5,-4,1}
};
vector<vector<int> > choice{
	{1},{1,3},{1,3},{1,3},{1,3},{1,3},{1,3,5},{1,2,3,4,5},{1,2,3,4,5,6}
};
vector<string> Roman{ "I","II","III","IV","V","VI" };
template<typename T>
vector<double> operator*(vector<T> v, double u)
{
	vector<double> ans;
	for (auto val : v)
		ans.push_back(val * u);
	return ans;
}
static void TranStrainToMatrix(vector<double> strain, vector<vector<double> >& matrix)
{
	matrix.resize(3, vector<double>(3));
	matrix[0][0] = strain[0] + 1;
	matrix[1][1] = strain[1] + 1;
	matrix[2][2] = strain[2] + 1;
	matrix[1][2] = matrix[2][1] = strain[3] / 2;
	matrix[0][2] = matrix[2][0] = strain[4] / 2;
	matrix[0][1] = matrix[1][0] = strain[5] / 2;
}
vector<string> vasp_file{ "INCAR","KPOINTS","POTCAR" };
elastic::elastic(vector<double> strain, const char file[])
{
	this->strain = strain;
	cal_method = { &elastic::cal_cub,&elastic::cal_hexa,&elastic::cal_trig1,&elastic::cal_trig2,
		&elastic::cal_tetra1,&elastic::cal_tetra2,&elastic::cal_ortho,&elastic::cal_mono,&elastic::cal_tric };
	Celas.resize(6, vector<double>(6));
	FILE* fp = fopen(file, "r");
	if (fp == nullptr)
	{
		cout << file << " IS NOT EXIST!" << endl;
		return;
	}
	readposcar(fp, pos);
	fclose(fp);
	double ntemp_vec[3][3];
	transpose_matrix(pos.vec, ntemp_vec);
	int* types = (int*)malloc(sizeof(int) * pos.nant[0]);
	translate_typenum_type(types, pos.nant[1], pos.typenum);
	SpglibDataset* dataset = spg_get_dataset(ntemp_vec, pos.xyz, types, pos.nant[0], SYMPREAC);
	int spacegroup_number = dataset->spacegroup_number;
	if (spacegroup_number == 1 || spacegroup_number == 2)
		Laue = N;
	if (spacegroup_number >= 3 && spacegroup_number < 16)
		Laue = M;
	if (spacegroup_number >= 16 && spacegroup_number < 75)
		Laue = O;
	if (spacegroup_number >= 75 && spacegroup_number < 89)
		Laue = T2;
	if (spacegroup_number >= 89 && spacegroup_number < 143)
		Laue = T1;
	if (spacegroup_number >= 143 && spacegroup_number < 149)
		Laue = R2;
	if (spacegroup_number >= 149 && spacegroup_number < 168)
		Laue = R1;
	if (spacegroup_number >= 168 && spacegroup_number < 195)
		Laue = H;
	if (spacegroup_number >= 195 && spacegroup_number < 230)
		Laue = C;
	lists = choice[Laue];
	stress.resize(lists.size(), vector<vector<double> >(strain.size(), vector<double>(6)));
	coeff.resize(lists.size(), vector<double>(6));
	spg_free_dataset(dataset);
	free(types);
}
void elastic::generatefile(const char filename[], vector<vector<double> > matrix)
{
	double n_vec[3][3];
	for (int i = 0; i < 3; i++)
	{
		double p[3] = { pos.vec[i][0],pos.vec[i][1],pos.vec[i][2] };
		for (int j = 0; j < 3; j++)
			n_vec[i][j] = p[0] * matrix[j][0] + p[1] * matrix[j][1] + p[2] * matrix[j][2];
	}
	POSCAR p(pos);
	for (int i = 0; i < 3; i++)
		for (int j = 0; j < 3; j++)
			p.vec[i][j] = n_vec[i][j];
	FILE* fp = fopen(filename, "w");
	savposcar(fp, p);
	fclose(fp);
}
void elastic::getstress(const char file[])
{
	vector<vector<vector<double> > > tmp(lists.size(), vector<vector<double> >(6, vector<double>(strain.size())));
	for (int i = 0; i < lists.size(); i++)
	{
		chdir(Roman[i].c_str());
		FILE* fp = fopen(file, "w");
		fprintf(fp, "	Strain       XX(GPa)       YY(GPa)       ZZ(GPa)       XY(GPa)       YZ(GPa)       ZX(GPa)\n");
		for (int j = 0; j < strain.size(); j++)
		{
			string foldname = "strain_" + to_string(strain[j]);
			chdir(foldname.c_str());
			stress[i][j] = get_stress() * (-0.1);
			fprintf(fp, "	%lf       %lf       %lf       %lf       %lf       %lf       %lf\n", strain[j], stress[i][j][0], stress[i][j][1],
				stress[i][j][2], stress[i][j][5], stress[i][j][3], stress[i][j][4]);
			chdir("..");
		}
		fclose(fp);
		chdir("..");
	}
	for (int i = 0; i < lists.size(); i++)
	{
		for (int j = 0; j < 6; j++)
		{
			for (int k = 0; k < stress[i].size(); k++)
				tmp[i][j][k] = stress[i][k][j];
			Eigen::VectorXd Coefficient(FitterLeastSquareMethod(strain, tmp[i][j], 1));
			coeff[i][j] = Coefficient[1];
		}
	}
}
void elastic::generate()
{
	for (int i = 0; i < lists.size(); i++)
	{
		string mkdir = "mkdir " + Roman[i];
		if (!access(Roman[i].c_str(), 0))
			cout << Roman[i] << " is exist! && skip this step!" << endl;
		else
			system(mkdir.c_str());
		cout << "Created " << Roman[i] << " folder!" << endl;
		chdir(Roman[i].c_str());
		for (auto val : strain)
		{
			vector<vector<double> > matrix;
			string foldname = "strain_" + to_string(val);
			string mkdir1 = "mkdir " + foldname;
			if (!access(foldname.c_str(), 0))
				cout << foldname << " is exist! && skip this step!" << endl;
			else
				system(mkdir1.c_str());
			cout << "Created " << foldname << " folder!" << endl;
			chdir(foldname.c_str());
			TranStrainToMatrix(deformation_type[lists[i] - 1] * val, matrix);
			generatefile("POSCAR", matrix);
			for (string file : vasp_file)
			{
				string source_file = "../../" + file;
				string dest_file = "./" + file;
				if (!access(source_file.c_str(), 0))
					copy_file(dest_file.c_str(), source_file.c_str());
			}
			chdir("..");
		}
		chdir("..");
	}
}
void elastic::calculate()
{
	cal_method[Laue](this);
	for (int i = 0; i < 6; i++)
		for (int j = 0; j < 6; j++)
			if (i > j)
				Celas[i][j] = Celas[j][i];
	PrintElasticProperty();
}
void elastic::PrintCelas(FILE* fp)
{
	for (int i = 0; i < 6; i++)
	{
		for (int j = 0; j < 6; j++)
			fprintf(fp, "	%lf ", Celas[i][j]);
		fprintf(fp, "\n");
	}
}
void elastic::PrintElasticProperty(const char file[])
{
	FILE* fp = fopen(file, "w");
	PrintCelas(fp);
	fclose(fp);
	cout << "Create " << file << " !" << endl;
}
void elastic::cal_cub()
{
	//%%%%%			Cubic system
	//%%%%% (3 independent elastic constants)
	//%%%%% c11  c12  c12    0    0    0
	//%%%%% c12  c11  c12    0    0    0
	//%%%%% c12  c12  c11    0    0    0
	//%%%%% 0    0    0		c44   0    0
	//%%%%% 0    0    0		 0   c44   0
	//%%%%% 0    0    0		 0    0   c44
	Eigen::Matrix<double, 3, 3> A;
	//c11 c12 c44
	A << 1, 5, 0,
		2, 4, 0,
		0, 0, 4;
	Eigen::Matrix<double, 3, 1> B;
	B << coeff[0][0], coeff[0][1], coeff[0][3];
	Eigen::Matrix<double, 3, 1> C(A.inverse() * B);
	Celas[0][0] = C[0]; //c11
	Celas[0][1] = C[1]; //c12
	Celas[0][2] = Celas[0][1]; //c13 = c12

	Celas[1][1] = Celas[0][0]; //c22 = c11
	Celas[1][2] = Celas[0][1]; //c23 = c12

	Celas[2][2] = Celas[0][0]; //c33 = c11

	Celas[3][3] = Celas[4][4] = Celas[5][5] = C[2]; //c44 = c55 = c66
}
void elastic::cal_hexa()
{
	//%%%%% Hexagonal system
	//%%%%% (5 independent elastic constants)
	//%%%%% c11  c12  c13    0    0    0
	//%%%%% c12  c11  c13    0    0    0
	//%%%%% c13  c13  c33    0    0    0
	//%%%%% 0    0    0    c44    0    0
	//%%%%% 0    0    0      0  c44    0
	//%%%%% 0    0    0      0    0  c66 = (c11 - c12) / 2
	Eigen::Matrix<double, 5, 5> A;
	//c11 c12 c13 c33 c44
	A << 1, 2, 3, 0, 0,
		2, 1, 3, 0, 0,
		0, 0, 3, 3, 0,
		0, 0, 0, 0, 4,
		3, -5, -1, 0, 0;
	Eigen::Matrix<double, 5, 1> B;
	B << coeff[0][0], coeff[0][1], coeff[0][2], coeff[0][3], coeff[1][0];
	Eigen::Matrix<double, 5, 1> C(A.inverse() * B);
	Celas[0][0] = C[0]; //c11
	Celas[0][1] = C[1]; //c12
	Celas[0][2] = C[2]; //c13

	Celas[1][1] = Celas[0][0]; //c22 = c12
	Celas[1][2] = Celas[0][2]; //c23 = c13

	Celas[2][2] = C[3]; //c33

	Celas[3][3] = Celas[4][4] = C[4]; //c44 = c55

	Celas[5][5] = (Celas[0][0] - Celas[0][1]) / 2; //c66
}
void elastic::cal_trig1()
{
	//%%%%%		Rhombohedral I system
	//%%%%% (6 independent elastic constants)
	//%%%%% c11  c12  c13  c14    0    0
	//%%%%% c12  c11  c13 -c14    0    0
	//%%%%% c13  c13  c33    0    0    0
	//%%%%% c14 -c14  0    c44    0    0
	//%%%%% 0    0    0    0     c44  c14
	//%%%%% 0    0    0    0     c14  c66 = (c11 - c12) / 2
	Eigen::Matrix<double, 6, 6> A;
	//c11 c12 c13 c14 c33 c44
	A << 1, 2, 3, 4, 0, 0,
		2, 1, 3, -4, 0, 0,
		0, 0, 3, 0, 3, 0,
		0, 0, 0, -1, 0, 4,
		0, 0, 0, 6, 0, 5,
		3, -5, -1, 6, 0, 0;
	Eigen::Matrix<double, 6, 1> B;
	B << coeff[0][0], coeff[0][1], coeff[0][2], coeff[0][3], coeff[0][4], coeff[1][0];
	Eigen::Matrix<double, 6, 1> C(A.inverse() * B);
	Celas[0][0] = C[0]; //c11
	Celas[0][1] = C[1]; //c12
	Celas[0][2] = C[2]; //c13
	Celas[0][3] = C[3]; //c14

	Celas[1][1] = Celas[0][0]; //c22 = c11
	Celas[1][2] = Celas[0][2]; //c23 = c13
	Celas[1][3] = -Celas[0][3]; //c24 = -c14

	Celas[2][2] = C[4]; //c33

	Celas[3][3] = C[5]; //c44

	Celas[4][4] = Celas[3][3]; //c55 = c44
	Celas[4][5] = Celas[0][3]; //c56 = c14

	Celas[5][5] = (Celas[0][0] - Celas[0][1]) / 2; //c66
}
void elastic::cal_trig2()
{
	//%%%%%		Rhombohedral II system
	//%%%%% (7 independent elastic constants)
	//%%%%% c11  c12  c13  c14  c15    0
	//%%%%% c12  c11  c13 -c14 -c15    0
	//%%%%% c13  c13  c33    0    0    0
	//%%%%% c14 -c14    0  c44    0 -c15
	//%%%%% c15 -c15    0    0  c44  c14
	//%%%%% 0    0      0 -c15  c14  c66 = (c11 - c12) / 2
	Eigen::Matrix<double, 7, 7> A;
	//c11 c12 c13 c14 c15 c33 c44
	A << 1, 2, 3, 4, 5, 0, 0,
		2, 1, 3, -4, -5, 0, 0,
		0, 0, 3, 0, 0, 3, 0,
		0, 0, 0, -1, -6, 0, 4,
		0, 0, 0, 0, 6, -1, 5,
		3, -5, -1, 6, 2, 0, 0,
		-5, 3, -1, -6, -2, 0, 0;
	Eigen::Matrix<double, 7, 1> B;
	B << coeff[0][0], coeff[0][1], coeff[0][2], coeff[0][3], coeff[0][4], coeff[1][0], coeff[1][1];
	Eigen::Matrix<double, 7, 1> C(A.inverse() * B);
	Celas[0][0] = C[0]; //c11
	Celas[0][1] = C[1]; //c12
	Celas[0][2] = C[2]; //c13
	Celas[0][3] = C[3]; //c14
	Celas[0][4] = C[4]; //c15

	Celas[1][1] = Celas[0][0]; //c22 = c11
	Celas[1][2] = Celas[0][2]; //c23 = c13
	Celas[1][3] = -Celas[0][3]; //c24 = -c14
	Celas[1][4] = -Celas[0][4]; //c25 = -c15

	Celas[2][2] = C[4]; //c33

	Celas[3][3] = C[5]; //c44
	Celas[3][5] = -Celas[0][4]; //c46 = -c15

	Celas[4][4] = Celas[3][3]; //c55 = c44
	Celas[4][5] = Celas[0][3]; //c56 = c14

	Celas[5][5] = (Celas[0][0] - Celas[0][1]) / 2; //c66
}
void elastic::cal_tetra1()
{
	//%%%%%		Tetragonal I system
	//%%%%% (6 independent elastic constants)
	//%%%%% c11  c12  c13    0    0    0
	//%%%%% c12  c11  c13    0    0    0
	//%%%%% c13  c13  c33    0    0    0
	//%%%%% 0    0    0     c44   0    0
	//%%%%% 0    0    0      0   c44   0
	//%%%%% 0    0    0      0    0   c66
	Eigen::Matrix<double, 6, 6> A;
	A << 1, 2, 3, 0, 0, 0,
		2, 1, 3, 0, 0, 0,
		0, 0, 3, 3, 0, 0,
		0, 0, 0, 0, 4, 0,
		0, 0, 0, 0, 0, 6,
		3, -5, -1, 0, 0, 0;
	Eigen::Matrix<double, 6, 1> B;
	B << coeff[0][0], coeff[0][1], coeff[0][2], coeff[0][3], coeff[0][5], coeff[1][0];
	Eigen::Matrix<double, 6, 1> C(A.inverse() * B);
	Celas[0][0] = C[0]; //c11
	Celas[0][1] = C[1]; //c12
	Celas[0][2] = C[2]; //c13

	Celas[1][1] = Celas[0][0]; //c22 = c12
	Celas[1][2] = Celas[0][2]; //c23 = c13

	Celas[2][2] = C[3]; //c33

	Celas[3][3] = Celas[4][4] = C[4]; //c44 = c55

	Celas[5][5] = C[5]; //c66
}
void elastic::cal_tetra2()
{
	//%%%%%		Tetragonal II system
	//%%%%% (7 independent elastic constants)
	//%%%%% c11  c12  c13    0    0   c16
	//%%%%% c12  c11  c13    0    0  -c16
	//%%%%% c13  c13  c33    0    0    0
	//%%%%% 0    0    0     c44   0    0
	//%%%%% 0    0    0      0   c44   0
	//%%%%% c16 -c16  0      0    0   c66
	Eigen::Matrix<double, 7, 7> A;
	A << 1, 2, 3, 0, 6, 0, 0,
		2, 1, 3, 0, -6, 0, 0,
		0, 0, 3, 3, 0, 0, 0,
		0, 0, 0, 0, 0, 4, 0,
		0, 0, 0, 0, 0, 0, 6,
		3, -5, -1, 0, -4, 0, 0,
		-5, 3, -1, 0, 4, 0, 0;
	Eigen::Matrix<double, 7, 1> B;
	B << coeff[0][0], coeff[0][1], coeff[0][2], coeff[0][3], coeff[0][5], coeff[1][0], coeff[1][1];
	Eigen::Matrix<double, 7, 1> C(A.inverse() * B);
	Celas[0][0] = C[0]; //c11
	Celas[0][1] = C[1]; //c12
	Celas[0][2] = C[2]; //c13
	Celas[0][5] = C[4]; //c16

	Celas[1][1] = Celas[0][0]; //c22 = c12
	Celas[1][2] = Celas[0][2]; //c23 = c13
	Celas[1][5] = -Celas[0][5]; //c26 = -c16

	Celas[2][2] = C[3]; //c33

	Celas[3][3] = Celas[4][4] = C[5]; //c44 = c55

	Celas[5][5] = C[6]; //c66
}
void elastic::cal_ortho()
{
	//%%%%%		Orthorhombic system
	//%%%%% (9 independent elastic constants)
	//%%%%% c11  c12  c13    0    0    0
	//%%%%%	c12  c22  c23    0    0    0
	//%%%%%	c13  c23  c33    0    0    0
	//%%%%%	0    0    0		c44   0    0
	//%%%%%	0    0    0		 0  c55    0
	//%%%%%	0    0    0		 0    0  c66
	Eigen::Matrix<double, 6, 6> _strain;
	Eigen::Matrix<double, 6, 6> _stress;
	_strain.fill(0);
	_stress.fill(0);
	for (int i = 0; i < 6; i++)
		for (int j = 0; j < 3; j++)
		{
			_strain(i, j) = deformation_type[lists[j] - 1][i];
			_stress(i, j) = coeff[j][i];
		}
	_strain(3, 3) = _strain(4, 4) = _strain(5, 5) = 1;
	_stress(3, 3) = coeff[0][3] / _strain(3,0);
	_stress(4, 4) = coeff[0][4] / _strain(4,0);
	_stress(5, 5) = coeff[0][5] / _strain(5,0);
	auto ans = _stress * (_strain.inverse());
	for (int i = 0; i < 6; i++)
		for (int j = 0; j < 6; j++)
			Celas[i][j] = ans(i,j);
}
void elastic::cal_mono()
{
	//[e14 e16 * [c44 c46 = [s14 s16
	//e24 e26]	 c46 c66]    s24 s26]
	Eigen::Matrix<double, 2, 2> tmp1; 
	tmp1 << deformation_type[lists[0] - 1][3], deformation_type[lists[0] - 1][5],
		deformation_type[lists[1] - 1][3], deformation_type[lists[1] - 1][5];
	//c66 c46
	Eigen::Matrix<double, 2, 2> tmp2;
	tmp2 << coeff[0][3], coeff[0][5],
		coeff[1][3], coeff[1][5];
	auto tmp3 = tmp1.inverse() * tmp2;
	double c44 = tmp3(0, 0);
	double c66 = tmp3(1, 1);
	Eigen::Matrix<double, 6, 6> _strain;
	Eigen::Matrix<double, 6, 6> _stress;
	_strain.fill(0);
	_stress.fill(0);
	for (int i = 0; i < 6; i++)
		for (int j = 0; j < 5; j++)
		{
			_strain(i, j) = deformation_type[lists[j] - 1][i];
			_stress(i, j) = coeff[j][i];
		}
	_strain(5, 5) = 1;
	_stress(5, 3) = c44;
	_stress(5, 5) = c66;
	auto ans = _stress * (_strain.inverse());
	for (int i = 0; i < 6; i++)
		for (int j = 0; j < 6; j++)
			Celas[i][j] = ans(i, j);
}
void elastic::cal_tric()
{
	Eigen::Matrix<double, 6, 6> _strain;
	Eigen::Matrix<double, 6, 6> _stress;
	for (int i = 0; i < 6; i++)
		for (int j = 0; j < 6; j++)
		{
			_strain(i, j) = deformation_type[lists[j] - 1][i];
			_stress(i, j) = coeff[j][i];
		}
	auto ans = _stress * (_strain.inverse());
	for (int i = 0; i < 6; i++)
		for (int j = 0; j < 6; j++)
			Celas[i][j] = ans(i, j);
}