#include "sgm_esl.h"

/**
 * @brief      初始化观测点数组
 * 
 * 注意我们的模型空间为z轴向下的右手直角坐标系，因此y轴的正方向为东向，而x轴的正方向为北向。这与通常二维数据的x与y轴朝向相反。
 *
 * @param      para_char  观测数据文件与参数
 *
 * @return     返回执行情况。正确返回0，否则返回-1
 */
int SGM_ESL::InitObsPoints(char* para_char){
	bool if_interpolate = false;
	char filename[1024];
	double obs_xmin, obs_xmax, obs_ymin, obs_ymax, obs_dx, obs_dy;
	int obs_xnum, obs_ynum;
	obspoint temp_obs;
	string temp_str;

	//解析文件名+r<xmin>/<xmax>/<ymin>/<ymax>+i<dx>/dy
	if (7 != sscanf(para_char, "%[^+]+r%lf/%lf/%lf/%lf+i%lf/%lf",
					filename, &obs_ymin, &obs_ymax, &obs_xmin, &obs_xmax, &obs_dy, &obs_dx)){
		strcpy(filename,para_char);
	}
	else{
		clog << BOLDYELLOW << "warning ==> " << RESET << "observations are recognized as a regular grid." << endl;
		if_interpolate = true;
		obs_ynum = (obs_ymax - obs_ymin)/obs_dy + 1;
		obs_xnum = (obs_xmax - obs_xmin)/obs_dx + 1;
	}
	
	ifstream infile;
	if (open_infile(infile,filename)) return -1;

	obspointArray temp_obs_points;
	while(getline(infile,temp_str)){
		if (*(temp_str.begin()) == '#') continue;
		else{
			//按每行5个数据解析 初始化为含观测值与不确定度的观测点
			if (5 == sscanf(temp_str.c_str(),"%lf %lf %lf %lf %lf",
				&temp_obs.y,&temp_obs.x,&temp_obs.z,&temp_obs.val,&temp_obs.dev)){
				temp_obs.z *= -1.0;
				temp_obs.id = temp_obs_points.size();
				temp_obs_points.push_back(temp_obs);
			}
			else{
				cerr << BOLDYELLOW << "ignored ==> " << RESET << "wrong input: " << temp_str << endl;
				continue;
			}
		}
	}
	infile.close();

	// 如果读入的观测点数为空则返回错误
	if (temp_obs_points.empty())
	{
		cerr << BOLDRED << "error ==> " << RESET << "fail to initial observations with the parameter: " << filename << endl;
		return -1;
	}
	// 根据所使用的等效源设置将观测数据插值至每个等效源块体的平面中心位置 以达到适当加密或抽稀观测数据的目的
	// 注意这里我们默认输入的观测点数据为规则网数据 数据的排序应为由左下到右上
	else if (if_interpolate){
		int temp_m, temp_n;
		obs_points_.reserve(mod_blocks_num_);
		for (int i = 0; i < mod_blocks_num_; i++){
			temp_obs.id = obs_points_.size();
			temp_obs.y = 0.5 * (mod_blocks_[i].ymin + mod_blocks_[i].ymax);
			temp_obs.x = 0.5 * (mod_blocks_[i].xmin + mod_blocks_[i].xmax);

			temp_m = floor((temp_obs.y - obs_ymin)/obs_dy);
			temp_n = floor((temp_obs.x - obs_xmin)/obs_dx);

			if (temp_m >= 0 && temp_n >= 0 && temp_m <= obs_ynum - 2 && temp_n <= obs_xnum - 2){
				temp_obs.z = cube_interpolate(obs_ymin + temp_m * obs_dy, obs_ymin + (temp_m + 1) * obs_dy,
											  obs_xmin + temp_n * obs_dx, obs_xmin + (temp_n + 1) * obs_dx,
											  temp_obs.y, temp_obs.x,
											  temp_obs_points[temp_n * obs_ynum + temp_m].z,
											  temp_obs_points[temp_n * obs_ynum + temp_m + 1].z,
											  temp_obs_points[(temp_n + 1) * obs_ynum + temp_m].z,
											  temp_obs_points[(temp_n + 1) * obs_ynum + temp_m + 1].z);

				temp_obs.val = cube_interpolate(obs_ymin + temp_m * obs_dy, obs_ymin + (temp_m + 1) * obs_dy,
											  obs_xmin + temp_n * obs_dx, obs_xmin + (temp_n + 1) * obs_dx,
											  temp_obs.y, temp_obs.x,
											  temp_obs_points[temp_n * obs_ynum + temp_m].val,
											  temp_obs_points[temp_n * obs_ynum + temp_m + 1].val,
											  temp_obs_points[(temp_n + 1) * obs_ynum + temp_m].val,
											  temp_obs_points[(temp_n + 1) * obs_ynum + temp_m + 1].val);

				temp_obs.dev = cube_interpolate(obs_ymin + temp_m * obs_dy, obs_ymin + (temp_m + 1) * obs_dy,
											  obs_xmin + temp_n * obs_dx, obs_xmin + (temp_n + 1) * obs_dx,
											  temp_obs.y, temp_obs.x,
											  temp_obs_points[temp_n * obs_ynum + temp_m].dev,
											  temp_obs_points[temp_n * obs_ynum + temp_m + 1].dev,
											  temp_obs_points[(temp_n + 1) * obs_ynum + temp_m].dev,
											  temp_obs_points[(temp_n + 1) * obs_ynum + temp_m + 1].dev);

				obs_points_.push_back(temp_obs);
			}
		}

		obs_points_num_ = obs_points_.size();
		//赋值给total_field_ 同时初始化其他数组大小
		total_field_.resize(obs_points_num_, 0.0);
		regional_field_.resize(obs_points_num_, 0.0);
		object_field_.resize(obs_points_num_, 0.0);
		regional_regional_field_.resize(obs_points_num_, 0.0);
		regional_object_field_.resize(obs_points_num_, 0.0);
		object_regional_field_.resize(obs_points_num_, 0.0);
		object_object_field_.resize(obs_points_num_, 0.0);
		for (int i = 0; i < obs_points_num_; i++){
			total_field_[i] = obs_points_[i].val;
		}
	}
	else{
		obs_points_ = temp_obs_points;

		obs_points_num_ = obs_points_.size();
		//赋值给total_field_ 同时初始化其他数组大小
		total_field_.resize(obs_points_num_,0.0);
		regional_field_.resize(obs_points_num_,0.0);
		object_field_.resize(obs_points_num_,0.0);
		regional_regional_field_.resize(obs_points_num_,0.0);
		regional_object_field_.resize(obs_points_num_,0.0);
		object_regional_field_.resize(obs_points_num_,0.0);
		object_object_field_.resize(obs_points_num_,0.0);
		for (int i = 0; i < obs_points_num_; i++){
			total_field_[i] = obs_points_[i].val;
		}
	}

	// 清理临时数组
	temp_obs_points.clear();
	vector<obspoint>().swap(temp_obs_points);
	return 0;
}