#include "xyz2shc.h"

void XYZ2SHC::CalCKernel(){
	int i,k,n;
	int start_index,former_index,flow_index;
	//初始化c_kernel_;
	c_kernel_.resize(obs_num_);
	for (int i = 0; i < obs_num_; i++){
		c_kernel_[i].resize(shc_num_,0.0);
	}

	//先计算对角线上的伴随勒让德系数 对角线上第一个值
#pragma omp parallel for private(i) schedule(guided)
	for (int i = 0; i < obs_num_; i++){
		c_kernel_[i][0] = 1.0;
		c_kernel_[i][half_shc_num_] = 1.0;
	}

	//对角线上第二个值
	start_index = shc_order_ + 1;
#pragma omp parallel for private(i) schedule(guided)
	for (int i = 0; i < obs_num_; i++){
		c_kernel_[i][start_index] = sqrt(3.0)*sin(obs_point_[i].theta*Pi/180.0);
		c_kernel_[i][start_index + half_shc_num_] = sqrt(3.0)*sin(obs_point_[i].theta*Pi/180.0);
	}

	//对角线上的其他值
	for (k = 2; k < shc_order_+1; k++){
		former_index = int(0.5*(k-1)*(2*shc_order_+4-k));
		start_index = int(0.5*k*(2*shc_order_+3-k));
#pragma omp parallel for private(i) shared(k,start_index,former_index) schedule(guided)
		for (int i = 0; i < obs_num_; i++){
			c_kernel_[i][start_index] = sin(obs_point_[i].theta*Pi/180.0)*sqrt(0.5*(2.0*k+1)/k)*c_kernel_[i][former_index];
			c_kernel_[i][start_index + half_shc_num_] = sin(obs_point_[i].theta*Pi/180.0)*sqrt(0.5*(2.0*k+1)/k)*c_kernel_[i][former_index + half_shc_num_];
		}
	}

	//计算副对角线上的值
	for (k = 0; k < shc_order_; k++){
		start_index = int(0.5*k*(2*shc_order_+3-k));
#pragma omp parallel for private(i) shared(k,start_index) schedule(guided)
		for (int i = 0; i < obs_num_; i++){
			c_kernel_[i][start_index + 1] = cos(obs_point_[i].theta*Pi/180.0)*sqrt(2.0*k+3)*c_kernel_[i][start_index];
			c_kernel_[i][start_index + half_shc_num_ + 1] = cos(obs_point_[i].theta*Pi/180.0)*sqrt(2.0*k+3)*c_kernel_[i][start_index + half_shc_num_];
		}
	}

	//计算其他位置上的值
	for (n = 0; n < shc_order_; n++){
		start_index = int(0.5*n*(2*shc_order_+3-n));
		for (k = n+2; k < shc_order_+1; k++){
			former_index = k - n;
#pragma omp parallel for private(i) shared(start_index,former_index,n,k) schedule(guided)
			for (int i = 0; i < obs_num_; i++){
				c_kernel_[i][start_index + former_index] = Anm_[k][n]*cos(obs_point_[i].theta*Pi/180.0)*c_kernel_[i][start_index+flow_index-1] - Bnm_[k][n]*c_kernel_[i][start_index+flow_index-2];
				c_kernel_[i][start_index + former_index + half_shc_num_] = Anm_[k][n]*cos(obs_point_[i].theta*Pi/180.0)*c_kernel_[i][start_index+flow_index-1+half_shc_num_] - Bnm_[k][n]*c_kernel_[i][start_index+flow_index-2+half_shc_num_];
			}
		}
	}

	//添加经度位置信息
	ProgressBar *bar = new ProgressBar(shc_order_+1,"Calculating Ckernel");
	for (n = 0; n < shc_order_+1; n++){
		bar->Progressed(i);

		start_index = int(0.5*n*(2*shc_order_+3-n));
		for (k = n; k < shc_order_+1; k++){
			former_index = k - n;
#pragma omp parallel for private(i) shared(start_index,former_index,n) schedule(guided)
			for (int i = 0; i < obs_num_; i++){
				c_kernel_[i][start_index + former_index] *= cos(n*obs_point_[i].phi*Pi/180.0);;
				c_kernel_[i][start_index + former_index + half_shc_num_] *= sin(n*obs_point_[i].phi*Pi/180.0);
			}
		}
	}
	return;
}