#include "XRD3D.cuh"

__global__ void Copy_Ucrd_From_System_To_XRD3D
(const int atom_numbers, const int *atom_serial, const UNSIGNED_INT_VECTOR *system_ucrd, UNSIGNED_INT_VECTOR *ucrd)
{
	int gpu_grid_layer_numbers = gridDim.x*blockDim.x;
	for (int i = blockIdx.x*blockDim.x + threadIdx.x; i < atom_numbers; i = i + gpu_grid_layer_numbers)
	{
		ucrd[i] = system_ucrd[atom_serial[i]];
	}
}
void XRD3D::Initial(CONTROLLER *controller,const int system_atom_numbers,const VECTOR boxlength, char *module_name)
{
	if (module_name == NULL)
	{
		strcpy(this->module_name, "XRD3D");
	}
	else
	{
		strcpy(this->module_name, module_name);
	}

	if (controller[0].Command_Exist(this->module_name,"atom_in_file"))
	{
		is_initialized = 0;
		controller[0].printf("START INITIALIZING XRD3D:\n");

		this->system_atom_numbers = system_atom_numbers;
		this->box_length = boxlength;
		volume = boxlength.x*boxlength.y*boxlength.z;

		//初始化格子尺度信息
		grid_dimension.int_x = Get_Fft_Patameter(2.*boxlength.x);
		grid_dimension.int_y = Get_Fft_Patameter(2.*boxlength.y);
		grid_dimension.int_z = Get_Fft_Patameter(2.*boxlength.z);
		if (controller[0].Command_Exist(this->module_name,"Nx"))
		{
			grid_dimension.int_x = atoi(controller[0].Command(this->module_name, "Nx"));
		}
		if (controller[0].Command_Exist(this->module_name, "Ny"))
		{
			grid_dimension.int_y = atoi(controller[0].Command(this->module_name, "Ny"));
		}
		if (controller[0].Command_Exist(this->module_name, "Nz"))
		{
			grid_dimension.int_z = atoi(controller[0].Command(this->module_name, "Nz"));
		}
		layer_numbers = grid_dimension.int_x*grid_dimension.int_y;
		grid_numbers = layer_numbers*grid_dimension.int_z;
		controller[0].printf("    box length:%f %f %f\ngrid numbers:%d %d %d = %d\n",
			boxlength.x,boxlength.y,boxlength.z,
			grid_dimension.int_x,grid_dimension.int_y,grid_dimension.int_z,
			grid_numbers);


		//读入xatom序号
		FILE *temp_in = NULL;
		Open_File_Safely(&temp_in, controller[0].Command(this->module_name, "atom_in_file"), "r");
		fscanf(temp_in, "%d", &xatom_numbers);
		if (xatom_numbers > system_atom_numbers)
		{
			controller[0].printf("	Error:XRD3D atom numbers %d is bigger than system atom numbers %d! Please Check!\n", xatom_numbers, system_atom_numbers);
			getchar();
		}
		Malloc_Safely((void**)&h_atom_serial, sizeof(int)*xatom_numbers);
		Cuda_Malloc_Safely((void**)&d_atom_serial, sizeof(int)*xatom_numbers);
		Malloc_Safely((void**)&h_fi, sizeof(float)*xatom_numbers);
		Cuda_Malloc_Safely((void**)&d_fi, sizeof(float)*xatom_numbers);
		Malloc_Safely((void**)&h_satom2xatom_serial, sizeof(int)*system_atom_numbers);
		Cuda_Malloc_Safely((void**)&d_satom2xatom_serial, sizeof(int)*system_atom_numbers);
		//体系全原子向需要xrd3d的原子的映射表
		for (int i = 0; i < system_atom_numbers; i = i + 1)
		{
			h_satom2xatom_serial[i] = -1;
		}
		for (int i = 0; i < xatom_numbers; i = i + 1)
		{
			fscanf(temp_in, "%d %f\n", &h_atom_serial[i], &h_fi[i]);
			h_satom2xatom_serial[h_atom_serial[i]] = i;
		}
		cudaMemcpy(d_satom2xatom_serial, h_satom2xatom_serial, sizeof(int)*system_atom_numbers, cudaMemcpyHostToDevice);
		cudaMemcpy(d_fi, h_fi, sizeof(float)*xatom_numbers, cudaMemcpyHostToDevice);
		cudaMemcpy(d_atom_serial, h_atom_serial, sizeof(int)*xatom_numbers, cudaMemcpyHostToDevice);
		
		Cuda_Malloc_Safely((void**)&atom_uint_crd, sizeof(UNSIGNED_INT_VECTOR)*xatom_numbers);
		Cuda_Malloc_Safely((void**)&atom_frc, sizeof(VECTOR)*xatom_numbers);
		Cuda_Malloc_Safely((void**)&system_atom_crd_grad, sizeof(VECTOR)*system_atom_numbers);
		Reset_List << <ceilf(3.*xatom_numbers / 128), 128 >> >(3 * xatom_numbers, (float*)atom_frc, 0.);
		Reset_List << <ceilf(3.*system_atom_numbers / 128), 128 >> >(3 * system_atom_numbers, (float*)system_atom_crd_grad, 0.);


		//获取xrd3d衍射光的初始化信息系
		if (controller[0].Command_Exist(this->module_name, "theta"))
		{
			theta = atof(controller[0].Command(this->module_name, "theta"));
		}
		if (controller[0].Command_Exist(this->module_name, "lambda"))
		{
			lambda = atof(controller[0].Command(this->module_name, "lambda"));
		}
		if (controller[0].Command_Exist(this->module_name, "attenuation_coefficient"))
		{
			s = atof(controller[0].Command(this->module_name, "attenuation_coefficient"));
		}
		k0 = Get_k(theta, lambda);
		controller[0].printf("    theta %f, lambda %f, attenuation_coefficient %f,\n",
			theta, lambda,s);
		controller[0].printf("    k0 %f\n",
			k0);

		//初始化计算倒空间求和所需要的内容
		fft.Initial(controller, this);

		//初始化能量相关
		Cuda_Malloc_Safely((void**)&d_sum_of_reciprocal_energy, sizeof(float));
		if (is_controller_printf_initialized == 0)
		{
			controller[0].Step_Print_Initial(this->module_name, "%.2f");
			is_controller_printf_initialized = 1;
		}
		controller[0].printf("    structure last modify date is %d\n", last_modify_date);
		controller[0].printf("END INITIALIZING XRD3D\n\n");
		is_initialized = 1;
	}
	else
	{

		is_initialized = 0;
		controller[0].printf("XRD3D IS NOT INITIALIZED\n\n");
	}
}
void XRD3D::Refresh_Volume(const VECTOR boxlength)
{
	this->box_length = boxlength;
	volume = boxlength.x*boxlength.y*boxlength.z;
	fft.Refresh_FFT_Volume();
}
void XRD3D::Calculate_XRD3D_And_Grad(const UNSIGNED_INT_VECTOR *uint_crd)
{
	if (is_initialized)
	{
		Copy_Ucrd_From_System_To_XRD3D << <20, 32 >> >(xatom_numbers, d_atom_serial, uint_crd, atom_uint_crd);
		fft.FFT_One_Time();
	}
}
float XRD3D::XRD3D_Cv(const UNSIGNED_INT_VECTOR *uint_crd,int is_download)
{
	if (is_initialized)
	{
		Copy_Ucrd_From_System_To_XRD3D << <20, 32 >> >(xatom_numbers, d_atom_serial, uint_crd, atom_uint_crd);
		fft.FFT_One_Time();
		if (is_download)
		{
			cudaMemcpy(&h_sum_of_reciprocal_energy, d_sum_of_reciprocal_energy, sizeof(float), cudaMemcpyDeviceToHost);
			return h_sum_of_reciprocal_energy;
		}
		else
		{
			return 0.;
		}
	}
	else
	{
		return 0.;
	}
}
float XRD3D::Get_k(float theta, float lambda)
{
	return 4 * CONSTANT_Pi / lambda * sinf(theta * CONSTANT_Pi / 180.);
}