#include "sgm_esl.h"
#include "progressBar.h"

void SGM_ESL::CalDeltaTKernel(){
	int i,j;
	//double k1,k2,k3,k4,k5,k6;
	double x1,x2,y1,y2,z1,z2;
	double w111,w112,w121,w122,w211,w212,w221,w222;
	double v111,v112,v121,v122,v211,v212,v221,v222;
	double R222,R122,R212,R112,R221,R121,R211,R111;
	double G222,G122,G212,G112,G221,G121,G211,G111;

	double *k1, *k2, *k3, *k4, *k5, *k6;
	k1 = new double [mod_blocks_num_];
	k2 = new double [mod_blocks_num_];
	k3 = new double [mod_blocks_num_];
	k4 = new double [mod_blocks_num_];
	k5 = new double [mod_blocks_num_];
	k6 = new double [mod_blocks_num_];

	double I, D;
	I0_ = I0_ * Pi / 180.0;
	D0_ = D0_ * Pi / 180.0;
	// 如果区域属性被激活 则使用区域磁化参数 否则使用局部磁化参数
	for (int i = 0; i < mod_blocks_num_; i++)
	{
		if (regional_model_list_[i] == 1.0)
		{
			I = reg_I_ * Pi / 180.0;
			D = reg_D_ * Pi / 180.0;
		}
		else
		{
			I = loc_I_ * Pi / 180.0;
			D = loc_D_ * Pi / 180.0;
		}

		k1[i]=cos(I0_)*sin(D0_)*sin(I)+sin(I0_)*cos(I)*sin(D);
		k2[i]=cos(I0_)*cos(D0_)*sin(I)+sin(I0_)*cos(I)*cos(D);
		k3[i]=cos(I0_)*cos(D0_)*cos(I)*sin(D)+cos(I0_)*sin(D0_)*cos(I)*cos(D);
		k4[i]=cos(I0_)*cos(D0_)*cos(I)*cos(D);
		k5[i]=cos(I0_)*sin(D0_)*cos(I)*sin(D);
		k6[i]=-sin(I0_)*sin(I);
	}

	if (GM_kernel_.empty())
	{
		GM_kernel_.resize(obs_points_num_);
		for (i = 0; i < obs_points_num_; i++){
			GM_kernel_[i].resize(mod_blocks_num_,0.0);
		}
	}

	ProgressBar *bar = new ProgressBar(obs_points_num_,"Calculating GMKernel");
	for (i = 0; i < obs_points_num_; i++){
		bar->Progressed(i);
#pragma omp parallel for private(j,R222,R122,R212,R112,R221,R121,R211,R111,\
		G222,G122,G212,G112,G221,G121,G211,G111,\
		w111,w112,w121,w122,w211,w212,w221,w222,\
		v111,v112,v121,v122,v211,v212,v221,v222,\
		x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
		for (j = 0; j < mod_blocks_num_; j++){
			x1 = mod_blocks_[j].xmin; x2 = mod_blocks_[j].xmax;
			y1 = mod_blocks_[j].ymin; y2 = mod_blocks_[j].ymax;
			z1 = mod_blocks_[j].zmin; z2 = mod_blocks_[j].zmax;

			R222=sqrt((x2-obs_points_[i].x)*(x2-obs_points_[i].x)+(y2-obs_points_[i].y)*(y2-obs_points_[i].y)+(z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			R122=sqrt((x1-obs_points_[i].x)*(x1-obs_points_[i].x)+(y2-obs_points_[i].y)*(y2-obs_points_[i].y)+(z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			R212=sqrt((x2-obs_points_[i].x)*(x2-obs_points_[i].x)+(y1-obs_points_[i].y)*(y1-obs_points_[i].y)+(z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			R112=sqrt((x1-obs_points_[i].x)*(x1-obs_points_[i].x)+(y1-obs_points_[i].y)*(y1-obs_points_[i].y)+(z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			R221=sqrt((x2-obs_points_[i].x)*(x2-obs_points_[i].x)+(y2-obs_points_[i].y)*(y2-obs_points_[i].y)+(z1-obs_points_[i].z)*(z1-obs_points_[i].z));
			R121=sqrt((x1-obs_points_[i].x)*(x1-obs_points_[i].x)+(y2-obs_points_[i].y)*(y2-obs_points_[i].y)+(z1-obs_points_[i].z)*(z1-obs_points_[i].z));
			R211=sqrt((x2-obs_points_[i].x)*(x2-obs_points_[i].x)+(y1-obs_points_[i].y)*(y1-obs_points_[i].y)+(z1-obs_points_[i].z)*(z1-obs_points_[i].z));
			R111=sqrt((x1-obs_points_[i].x)*(x1-obs_points_[i].x)+(y1-obs_points_[i].y)*(y1-obs_points_[i].y)+(z1-obs_points_[i].z)*(z1-obs_points_[i].z));

			w222=((x2-obs_points_[i].x)*(x2-obs_points_[i].x))+R222*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			w122=((x1-obs_points_[i].x)*(x1-obs_points_[i].x))+R122*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			w212=((x2-obs_points_[i].x)*(x2-obs_points_[i].x))+R212*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			w112=((x1-obs_points_[i].x)*(x1-obs_points_[i].x))+R112*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			w221=((x2-obs_points_[i].x)*(x2-obs_points_[i].x))+R221*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
			w121=((x1-obs_points_[i].x)*(x1-obs_points_[i].x))+R121*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
			w211=((x2-obs_points_[i].x)*(x2-obs_points_[i].x))+R211*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
			w111=((x1-obs_points_[i].x)*(x1-obs_points_[i].x))+R111*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));

			v222=((y2-obs_points_[i].y)*(y2-obs_points_[i].y))+R222*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			v122=((y2-obs_points_[i].y)*(y2-obs_points_[i].y))+R122*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			v212=((y1-obs_points_[i].y)*(y1-obs_points_[i].y))+R212*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			v112=((y1-obs_points_[i].y)*(y1-obs_points_[i].y))+R112*(z2-obs_points_[i].z)+((z2-obs_points_[i].z)*(z2-obs_points_[i].z));
			v221=((y2-obs_points_[i].y)*(y2-obs_points_[i].y))+R221*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
			v121=((y2-obs_points_[i].y)*(y2-obs_points_[i].y))+R121*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
			v211=((y1-obs_points_[i].y)*(y1-obs_points_[i].y))+R211*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));
			v111=((y1-obs_points_[i].y)*(y1-obs_points_[i].y))+R111*(z1-obs_points_[i].z)+((z1-obs_points_[i].z)*(z1-obs_points_[i].z));

			G222=k1[j]*log(R222+x2-obs_points_[i].x)+k2[j]*log(R222+y2-obs_points_[i].y)+k3[j]*log(R222+z2-obs_points_[i].z)
			+k4[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/w222)+k5[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/v222)
			+k6[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/(R222*(z2-obs_points_[i].z)));

			G122=k1[j]*log(R122+x1-obs_points_[i].x)+k2[j]*log(R122+y2-obs_points_[i].y)+k3[j]*log(R122+z2-obs_points_[i].z)
			+k4[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/w122)+k5[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/v122)
			+k6[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/(R122*(z2-obs_points_[i].z)));

			G212=k1[j]*log(R212+x2-obs_points_[i].x)+k2[j]*log(R212+y1-obs_points_[i].y)+k3[j]*log(R212+z2-obs_points_[i].z)
			+k4[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/w212)+k5[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/v212)
			+k6[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/(R212*(z2-obs_points_[i].z)));

			G112=k1[j]*log(R112+x1-obs_points_[i].x)+k2[j]*log(R112+y1-obs_points_[i].y)+k3[j]*log(R112+z2-obs_points_[i].z)
			+k4[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/w112)+k5[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/v112)
			+k6[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/(R112*(z2-obs_points_[i].z)));

			G221=k1[j]*log(R221+x2-obs_points_[i].x)+k2[j]*log(R221+y2-obs_points_[i].y)+k3[j]*log(R221+z1-obs_points_[i].z)
			+k4[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/w221)+k5[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/v221)
			+k6[j]*atan(((x2-obs_points_[i].x)*(y2-obs_points_[i].y))/(R221*(z1-obs_points_[i].z)));

			G121=k1[j]*log(R121+x1-obs_points_[i].x)+k2[j]*log(R121+y2-obs_points_[i].y)+k3[j]*log(R121+z1-obs_points_[i].z)
			+k4[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/w121)+k5[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/v121)
			+k6[j]*atan(((x1-obs_points_[i].x)*(y2-obs_points_[i].y))/(R121*(z1-obs_points_[i].z)));

			G211=k1[j]*log(R211+x2-obs_points_[i].x)+k2[j]*log(R211+y1-obs_points_[i].y)+k3[j]*log(R211+z1-obs_points_[i].z)
			+k4[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/w211)+k5[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/v211)
			+k6[j]*atan(((x2-obs_points_[i].x)*(y1-obs_points_[i].y))/(R211*(z1-obs_points_[i].z)));

			G111=k1[j]*log(R111+x1-obs_points_[i].x)+k2[j]*log(R111+y1-obs_points_[i].y)+k3[j]*log(R111+z1-obs_points_[i].z)
			+k4[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/w111)+k5[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/v111)
			+k6[j]*atan(((x1-obs_points_[i].x)*(y1-obs_points_[i].y))/(R111*(z1-obs_points_[i].z)));

			GM_kernel_[i][j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111);
		}
	}

	delete[] k1;
	delete[] k2;
	delete[] k3;
	delete[] k4;
	delete[] k5;
	delete[] k6;
	return;
}