initial upload

.gitignore

@@ -30,3 +30,6 @@
 *.exe
 *.out
 *.app
+
+build/
+.DS_Store

CMakeLists.txt

@@ -0,0 +1,3 @@
+cmake_minimum_required(VERSION 3.15.2)
+project(GM3D)
+add_subdirectory(src)

src/CMakeLists.txt

@@ -0,0 +1,10 @@
+aux_source_directory(. SRC_DIR)
+
+set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin)
+
+add_executable(gm3d ${SRC_DIR})
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --std=c++11 -O2")
+
+set(CMAKE_INSTALL_PREFIX /usr/local)
+install(TARGETS gm3d RUNTIME DESTINATION sbin)

src/add_interface_block.cpp

@@ -0,0 +1,145 @@
+#include "gm3d.h"
+//我们读入一个界面数据 插值计算每个块体中心位置的值 然后按情况赋值
+int GM3D::AddInterfaceBlock(modelist para_list){
+	int m, n, xnum, ynum;
+	double xs,xe,xmin,xmax,dx;
+	double ys,ye,ymin,ymax,dy;
+	bool model_added = false;
+
+	string temp_str;
+	stringstream temp_ss;
+
+	cpoint temp_topo;
+	cpointArray face_topo;
+
+	ifstream infile;
+	char filename[1024];
+	if (1 == sscanf(para_list.mod_para,"%s",filename)){
+		if (open_infile(infile,filename)) return -1;
+
+		while (getline(infile,temp_str)){
+			//#range必须出现在数据之前
+			if (*(temp_str.begin()) == '#'){
+				if (6 == sscanf(temp_str.c_str(),"# range=%lf/%lf/%lf/%lf/%lf/%lf",&xs,&dx,&xe,&ys,&dy,&ye)){
+					xmin = MIN(xs,xe); xmax = MAX(xs,xe);
+					ymin = MIN(ys,ye); ymax = MAX(ys,ye);
+					dx = fabs(dx); dy = fabs(dy);
+					xnum = round((xmax - xmin)/dx) + 1;
+					ynum = round((ymax - ymin)/dy) + 1;
+				}
+				else continue;
+			}
+			else{
+				temp_ss = str2ss(temp_str);
+				temp_ss >> temp_topo.x >> temp_topo.y >> temp_topo.z;
+				face_topo.push_back(temp_topo);
+			}
+		}
+		infile.close();
+
+		if (!strcmp(para_list.val_type,"replace/top")){
+			for (int i = 0; i < model_num_; i++){
+				m = floor((model_cube_[i].cen.x - xmin)/dx);
+				n = floor((model_cube_[i].cen.y - ymin)/dy);
+
+				temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
+											face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
+
+				if (model_cube_[i].cen.z >= temp_topo.z){
+					model_block_val_[i] = para_list.mod_value;
+					model_added = true;
+				}
+			}
+		}
+		else if (!strcmp(para_list.val_type,"replace/bot")){
+			for (int i = 0; i < model_num_; i++){
+				m = floor((model_cube_[i].cen.x - xmin)/dx);
+				n = floor((model_cube_[i].cen.y - ymin)/dy);
+
+				temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
+											face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
+
+				if (model_cube_[i].cen.z <= temp_topo.z){
+					model_block_val_[i] = para_list.mod_value;
+					model_added = true;
+				}
+			}
+		}
+		else if (!strcmp(para_list.val_type,"add/top")){
+			for (int i = 0; i < model_num_; i++){
+				m = floor((model_cube_[i].cen.x - xmin)/dx);
+				n = floor((model_cube_[i].cen.y - ymin)/dy);
+
+				temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
+											face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
+
+				if (model_cube_[i].cen.z >= temp_topo.z){
+					if (model_block_val_[i] == BDL_MAX)
+						model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
+					else
+						model_block_val_[i] += para_list.mod_value;
+					model_added = true;
+				}
+			}
+		}
+		else if (!strcmp(para_list.val_type,"add/bot")){
+			for (int i = 0; i < model_num_; i++){
+				m = floor((model_cube_[i].cen.x - xmin)/dx);
+				n = floor((model_cube_[i].cen.y - ymin)/dy);
+
+				temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
+											face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
+
+				if (model_cube_[i].cen.z <= temp_topo.z){
+					if (model_block_val_[i] == BDL_MAX)
+						model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
+					else
+						model_block_val_[i] += para_list.mod_value;
+					model_added = true;
+				}
+			}
+		}
+		else if (!strcmp(para_list.val_type,"erase/top")){
+			for (int i = 0; i < model_num_; i++){
+				m = floor((model_cube_[i].cen.x - xmin)/dx);
+				n = floor((model_cube_[i].cen.y - ymin)/dy);
+
+				temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
+											face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
+
+				if (model_cube_[i].cen.z >= temp_topo.z){
+					model_block_val_[i] = BDL_MAX;
+					model_added = true;
+				}
+			}
+		}
+		else if (!strcmp(para_list.val_type,"erase/bot")){
+			for (int i = 0; i < model_num_; i++){
+				m = floor((model_cube_[i].cen.x - xmin)/dx);
+				n = floor((model_cube_[i].cen.y - ymin)/dy);
+
+				temp_topo.z = grid_interpolate(xmin+dx*m,ymin+dy*n,dx,dy,model_cube_[i].cen.x,model_cube_[i].cen.y,
+											face_topo[xnum*n+m].z,face_topo[xnum*n+m+1].z,face_topo[xnum*(n+1)+m].z,face_topo[xnum*(n+1)+m+1].z);
+
+				if (model_cube_[i].cen.z <= temp_topo.z){
+					model_block_val_[i] = BDL_MAX;
+					model_added = true;
+				}
+			}
+		}
+		else{
+			cerr << BOLDRED << "error ==> " << RESET << "wrong value type: " << para_list.val_type << " of the model type: " << para_list.mod_type << endl;
+			return -1;
+		}
+	}
+	else{
+		cerr << BOLDRED << "error ==> " << RESET << "fail to add blocks with the parameter: " << para_list.mod_para << endl;
+		return -1;
+	}
+
+	if (!model_added){
+		cerr << BOLDYELLOW << "warning ==> " << RESET << "no block changed with the parameter: " << para_list.mod_para << endl;
+		return -1;
+	}
+	return 0;
+}

src/add_models.cpp

@@ -0,0 +1,45 @@
+#include "gm3d.h"
+
+int GM3D::AddModels(char* filename){
+	string temp_str;
+	stringstream temp_ss;
+	modelist temp_list;
+
+	ifstream infile;
+	if (open_infile(infile,filename)) return -1;
+	while(getline(infile,temp_str)){
+		if (*(temp_str.begin()) == '#') continue;
+		else{
+			//按每行5个数据解析 初始化为含观测值与不确定度的观测点
+			if (4 == sscanf(temp_str.c_str(),"%s %s %lf %s",
+				temp_list.mod_type,temp_list.val_type,&temp_list.mod_value,temp_list.mod_para)){
+				model_list_.push_back(temp_list);
+			}
+			else{
+				cerr << BOLDYELLOW << "ignored ==> " << RESET << "wrong input: " << temp_str << endl;
+				continue;
+			}
+		}
+	}
+	infile.close();
+
+	for (int i = 0; i < model_list_.size(); i++){
+		if (!strcmp(model_list_[i].mod_type,"regular_block")){
+			AddRegularBlock(model_list_[i]);
+		}
+		else if (!strcmp(model_list_[i].mod_type,"tilted_block")){
+			AddTiltedBlock(model_list_[i]);
+		}
+		else if (!strcmp(model_list_[i].mod_type,"sphere")){
+			AddSphereBlock(model_list_[i]);
+		}
+		else if (!strcmp(model_list_[i].mod_type,"interface")){
+			AddInterfaceBlock(model_list_[i]);
+		}
+		else{
+			cerr << BOLDYELLOW << "ignored ==> " << RESET << "unknown model type: " << model_list_[i].mod_type << endl;
+			continue;
+		}
+	}
+	return 0;
+}

src/add_regular_block.cpp

@@ -0,0 +1,62 @@
+#include "gm3d.h"
+
+int GM3D::AddRegularBlock(modelist para_list){
+	double xs,xe,xmin,xmax;
+	double ys,ye,ymin,ymax;
+	double zs,ze,zmin,zmax;
+	bool model_added = false;
+
+	if (6 == sscanf(para_list.mod_para,"%lf/%lf/%lf/%lf/%lf/%lf",&xs,&xe,&ys,&ye,&zs,&ze)){
+		xmin = MIN(xs,xe); xmax = MAX(xs,xe);
+		ymin = MIN(ys,ye); ymax = MAX(ys,ye);
+		zmin = MIN(zs,ze); zmax = MAX(zs,ze);
+
+		if (!strcmp(para_list.val_type,"replace")){
+			for (int i = 0; i < model_num_; i++){
+				if (model_cube_[i].cen.x >= xmin && model_cube_[i].cen.x <= xmax &&
+					model_cube_[i].cen.y >= ymin && model_cube_[i].cen.y <= ymax &&
+					model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
+					model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
+					model_added = true;
+				}
+			}
+		}
+		else if (!strcmp(para_list.val_type,"add")){
+			for (int i = 0; i < model_num_; i++){
+				if (model_cube_[i].cen.x >= xmin && model_cube_[i].cen.x <= xmax &&
+					model_cube_[i].cen.y >= ymin && model_cube_[i].cen.y <= ymax &&
+					model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
+					if (model_block_val_[i] == BDL_MAX)
+						model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
+					else
+						model_block_val_[i] += para_list.mod_value;
+					model_added = true;
+				}
+			}
+		}
+		else if (!strcmp(para_list.val_type,"erase")){
+			for (int i = 0; i < model_num_; i++){
+				if (model_cube_[i].cen.x >= xmin && model_cube_[i].cen.x <= xmax &&
+					model_cube_[i].cen.y >= ymin && model_cube_[i].cen.y <= ymax &&
+					model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
+					model_block_val_[i] = BDL_MAX; //注意重复赋值的块体会覆盖
+					model_added = true;
+				}
+			}
+		}
+		else{
+			cerr << BOLDRED << "error ==> " << RESET << "wrong value type: " << para_list.val_type << " of the model type: " << para_list.mod_type << endl;
+			return -1;
+		}
+	}
+	else{
+		cerr << BOLDRED << "error ==> " << RESET << "fail to add blocks with the parameter: " << para_list.mod_para << endl;
+		return -1;
+	}
+
+	if (!model_added){
+		cerr << BOLDYELLOW << "warning ==> " << RESET << "no block changed with the parameter: " << para_list.mod_para << endl;
+		return -1;
+	}
+	return 0;
+}

src/add_sphere_block.cpp

@@ -0,0 +1,82 @@
+#include "gm3d.h"
+
+int GM3D::AddSphereBlock(modelist para_list){
+	double xc,yc,zc,rad_x,rad_y,rad_z;
+	double dist,rad_limit;
+	double rel_x,rel_y,rel_z,theta,phi;
+	bool model_added = false;
+
+	if (6 == sscanf(para_list.mod_para,"%lf/%lf/%lf/%lf/%lf/%lf",&xc,&yc,&zc,&rad_x,&rad_y,&rad_z)){
+		if (!strcmp(para_list.val_type,"replace")){
+			for (int i = 0; i < model_num_; i++){
+				rel_x = model_cube_[i].cen.x - xc;
+				rel_y = model_cube_[i].cen.y - yc;
+				rel_z = model_cube_[i].cen.z - zc;
+				dist = sqrt(rel_x*rel_x + rel_y*rel_y + rel_z*rel_z);
+
+				theta = acos(rel_z/dist);
+				phi = atan2(rel_y,rel_x);
+
+				rad_limit = rad_x*rad_y*rad_z/sqrt(pow(rad_y*rad_z*sin(theta)*cos(phi),2) + pow(rad_x*rad_z*sin(theta)*sin(phi),2) + pow(rad_x*rad_y*cos(theta),2));
+
+				if (dist <= rad_limit){
+					model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
+					model_added = true;
+				}
+			}
+		}
+		else if (!strcmp(para_list.val_type,"add")){
+			for (int i = 0; i < model_num_; i++){
+				rel_x = model_cube_[i].cen.x - xc;
+				rel_y = model_cube_[i].cen.y - yc;
+				rel_z = model_cube_[i].cen.z - zc;
+				dist = sqrt(rel_x*rel_x + rel_y*rel_y + rel_z*rel_z);
+
+				theta = acos(rel_z/dist);
+				phi = atan2(rel_y,rel_x);
+
+				rad_limit = rad_x*rad_y*rad_z/sqrt(pow(rad_y*rad_z*sin(theta)*cos(phi),2) + pow(rad_x*rad_z*sin(theta)*sin(phi),2) + pow(rad_x*rad_y*cos(theta),2));
+
+				if (dist <= rad_limit){
+					if (model_block_val_[i] == BDL_MAX)
+						model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
+					else
+						model_block_val_[i] += para_list.mod_value;
+					model_added = true;
+				}
+			}
+		}
+		else if (!strcmp(para_list.val_type,"erase")){
+			for (int i = 0; i < model_num_; i++){
+				rel_x = model_cube_[i].cen.x - xc;
+				rel_y = model_cube_[i].cen.y - yc;
+				rel_z = model_cube_[i].cen.z - zc;
+				dist = sqrt(rel_x*rel_x + rel_y*rel_y + rel_z*rel_z);
+
+				theta = acos(rel_z/dist);
+				phi = atan2(rel_y,rel_x);
+
+				rad_limit = rad_x*rad_y*rad_z/sqrt(pow(rad_y*rad_z*sin(theta)*cos(phi),2) + pow(rad_x*rad_z*sin(theta)*sin(phi),2) + pow(rad_x*rad_y*cos(theta),2));
+
+				if (dist <= rad_limit){
+					model_block_val_[i] = BDL_MAX; //注意重复赋值的块体会覆盖
+					model_added = true;
+				}
+			}
+		}
+		else{
+			cerr << BOLDRED << "error ==> " << RESET << "wrong value type: " << para_list.val_type << " of the model type: " << para_list.mod_type << endl;
+			return -1;
+		}
+	}
+	else{
+		cerr << BOLDRED << "error ==> " << RESET << "fail to add blocks with the parameter: " << para_list.mod_para << endl;
+		return -1;
+	}
+
+	if (!model_added){
+		cerr << BOLDYELLOW << "warning ==> " << RESET << "no block changed with the parameter: " << para_list.mod_para << endl;
+		return -1;
+	}
+	return 0;
+}

src/add_tilted_block.cpp

@@ -0,0 +1,86 @@
+#include "gm3d.h"
+
+int GM3D::AddTiltedBlock(modelist para_list){
+	double xs_1,xe_1,xmin_1,xmax_1;
+	double xs_2,xe_2,xmin_2,xmax_2;
+	double ys_1,ye_1,ymin_1,ymax_1;
+	double ys_2,ye_2,ymin_2,ymax_2;
+	double layer_xmin,layer_xmax,layer_ymin,layer_ymax;
+	double zs,ze,zmin,zmax;
+	bool model_added = false;
+
+	if (10 == sscanf(para_list.mod_para,"%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf",
+		&xs_1,&xe_1,&ys_1,&ye_1,&zs,&xs_2,&xe_2,&ys_2,&ye_2,&ze)){
+		xmin_1 = MIN(xs_1,xe_1); xmax_1 = MAX(xs_1,xe_1);
+		ymin_1 = MIN(ys_1,ye_1); ymax_1 = MAX(ys_1,ye_1);
+		xmin_2 = MIN(xs_2,xe_2); xmax_2 = MAX(xs_2,xe_2);
+		ymin_2 = MIN(ys_2,ye_2); ymax_2 = MAX(ys_2,ye_2);
+		zmin = MIN(zs,ze); zmax = MAX(zs,ze);
+
+		if (!strcmp(para_list.val_type,"replace")){
+			for (int i = 0; i < model_num_; i++){
+				//计算当前层的x与y范围
+				layer_xmin = (model_cube_[i].cen.z - zmin)*(xmin_2 - xmin_1)/(zmax - zmin) + xmin_1;
+				layer_xmax = (model_cube_[i].cen.z - zmin)*(xmax_2 - xmax_1)/(zmax - zmin) + xmax_1;
+				layer_ymin = (model_cube_[i].cen.z - zmin)*(ymin_2 - ymin_1)/(zmax - zmin) + ymin_1;
+				layer_ymax = (model_cube_[i].cen.z - zmin)*(ymax_2 - ymax_1)/(zmax - zmin) + ymax_1;
+
+				if (model_cube_[i].cen.x >= layer_xmin && model_cube_[i].cen.x <= layer_xmax &&
+					model_cube_[i].cen.y >= layer_ymin && model_cube_[i].cen.y <= layer_ymax &&
+					model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
+					model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
+					model_added = true;
+				}
+			}
+		}
+		else if (!strcmp(para_list.val_type,"add")){
+			for (int i = 0; i < model_num_; i++){
+				//计算当前层的x与y范围
+				layer_xmin = (model_cube_[i].cen.z - zmin)*(xmin_2 - xmin_1)/(zmax - zmin) + xmin_1;
+				layer_xmax = (model_cube_[i].cen.z - zmin)*(xmax_2 - xmax_1)/(zmax - zmin) + xmax_1;
+				layer_ymin = (model_cube_[i].cen.z - zmin)*(ymin_2 - ymin_1)/(zmax - zmin) + ymin_1;
+				layer_ymax = (model_cube_[i].cen.z - zmin)*(ymax_2 - ymax_1)/(zmax - zmin) + ymax_1;
+
+				if (model_cube_[i].cen.x >= layer_xmin && model_cube_[i].cen.x <= layer_xmax &&
+					model_cube_[i].cen.y >= layer_ymin && model_cube_[i].cen.y <= layer_ymax &&
+					model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
+					if (model_block_val_[i] == BDL_MAX)
+						model_block_val_[i] = para_list.mod_value; //注意重复赋值的块体会覆盖
+					else
+						model_block_val_[i] += para_list.mod_value;
+					model_added = true;
+				}
+			}
+		}
+		else if (!strcmp(para_list.val_type,"erase")){
+			for (int i = 0; i < model_num_; i++){
+				//计算当前层的x与y范围
+				layer_xmin = (model_cube_[i].cen.z - zmin)*(xmin_2 - xmin_1)/(zmax - zmin) + xmin_1;
+				layer_xmax = (model_cube_[i].cen.z - zmin)*(xmax_2 - xmax_1)/(zmax - zmin) + xmax_1;
+				layer_ymin = (model_cube_[i].cen.z - zmin)*(ymin_2 - ymin_1)/(zmax - zmin) + ymin_1;
+				layer_ymax = (model_cube_[i].cen.z - zmin)*(ymax_2 - ymax_1)/(zmax - zmin) + ymax_1;
+
+				if (model_cube_[i].cen.x >= layer_xmin && model_cube_[i].cen.x <= layer_xmax &&
+					model_cube_[i].cen.y >= layer_ymin && model_cube_[i].cen.y <= layer_ymax &&
+					model_cube_[i].cen.z >= zmin && model_cube_[i].cen.z <= zmax){
+					model_block_val_[i] = BDL_MAX; //注意重复赋值的块体会覆盖
+					model_added = true;
+				}
+			}
+		}
+		else{
+			cerr << BOLDRED << "error ==> " << RESET << "wrong value type: " << para_list.val_type << " of the model type: " << para_list.mod_type << endl;
+			return -1;
+		}
+	}
+	else{
+		cerr << BOLDRED << "error ==> " << RESET << "fail to add blocks with the parameter: " << para_list.mod_para << endl;
+		return -1;
+	}
+
+	if (!model_added){
+		cerr << BOLDYELLOW << "warning ==> " << RESET << "no block changed with the parameter: " << para_list.mod_para << endl;
+		return -1;
+	}
+	return 0;
+}

src/build_regular_grid.cpp

@@ -0,0 +1,70 @@
+#include "gm3d.h"
+
+int GM3D::BuildRegularGrid(char* space_para){
+	cpoint temp_cp;
+	cube temp_cu;
+	string temp_id_str;
+	double x,xs,dx,xe,xmin,xmax;
+	double y,ys,dy,ye,ymin,ymax;
+	double z,zs,dz,ze,zmin,zmax;
+	double sign[8][3] = {{-0.5,-0.5,-0.5},{0.5,-0.5,-0.5},{0.5,0.5,-0.5},{-0.5,0.5,-0.5},
+						{-0.5,-0.5,0.5},{0.5,-0.5,0.5},{0.5,0.5,0.5},{-0.5,0.5,0.5}};
+
+	_str2pointMap map_str_point;
+	_str2pointMap::iterator imsp;
+
+	if (9 == sscanf(space_para,"%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf",
+					&xs,&dx,&xe,&ys,&dy,&ye,&zs,&dz,&ze)){
+		xmin = MIN(xs,xe); xmax = MAX(xs,xe);
+		ymin = MIN(ys,ye); ymax = MAX(ys,ye);
+		zmin = MIN(zs,ze); zmax = MAX(zs,ze);
+
+		temp_cu.dx = fabs(dx); temp_cu.dy = fabs(dy); temp_cu.dz = fabs(dz);
+		y = ys;
+		while(y >= ymin && y <= ymax){
+			x = xs;
+			while(x >= xmin && x <= xmax){
+				z = zs;
+				while(z >= zmin && z <= zmax){
+					//添加denMod
+					temp_cu.cen.id = model_cube_.size();
+					temp_cu.cen.x = x; temp_cu.cen.y = y; temp_cu.cen.z = z;
+					//添加mshVert
+					for (int i = 0; i < 8; i++){
+						temp_cp.id = model_vert_.size(); //添加msh的顶点索引为mshVert的大小
+						temp_cp.x = temp_cu.cen.x - sign[i][0]*temp_cu.dx; //左下底角
+						temp_cp.y = temp_cu.cen.y - sign[i][1]*temp_cu.dy;
+						temp_cp.z = temp_cu.cen.z - sign[i][2]*temp_cu.dz;
+						temp_id_str = cpoint_id(temp_cp);
+						imsp = map_str_point.find(temp_id_str);
+						//利用map_vert查到当前顶点是否存在,这里需要注意，如果顶点已经存在则只需要将顶点索引置为已存在顶点的索引，不增加顶点计数
+						if(imsp!=map_str_point.end()){
+							temp_cu.ids[i] = imsp->second.id;
+						}
+						//若为新的顶点则将其增加到两个映射和一个链表中
+						else{
+							temp_cu.ids[i] = temp_cp.id;//新的顶点索引等于顶点集的数量
+							model_vert_.push_back(temp_cp);//将新产生的顶点保存到顶点链表中
+							map_str_point[temp_id_str] = temp_cp;//将新产生的顶点保存到顶点位置映射中
+						}
+					}
+					model_cube_.push_back(temp_cu);
+					z += dz;
+				}
+				x += dx;
+			}
+			y += dy;
+		}
+	}
+
+	if (model_cube_.empty()){
+		cerr << BOLDRED << "error ==> " << RESET << "fail to initial model space with the parameter: " << space_para << endl;
+		return -1;
+	}
+	else{
+		vert_num_ = model_vert_.size();
+		model_num_ = model_cube_.size();
+		model_block_val_.resize(model_num_,BDL_MAX); //初始化模型块体值为BDL_MAX
+	}
+	return 0;
+}

src/disp_help.cpp

@@ -0,0 +1,329 @@
+#include "disp_help.h"
+
+void DispHelp::addHeadInfo(string s1,string s2,string s3,string s4)
+{
+	ex_name = s1; version = s2; descript = s3; author = s4;
+	return;
+}
+
+void DispHelp::addUsage(string usg)
+{
+	usages.push_back(usg);
+	return;
+}
+
+void DispHelp::addOption(string msg,string sflag,string lflag)
+{
+	option tmp_option;
+	tmp_option.message = msg; tmp_option.flag_s = sflag; tmp_option.flag_l = lflag;
+	options.push_back(tmp_option);
+	return;
+}
+
+void DispHelp::addOptionSec(string msg,int index)
+{
+	if (index < 0)
+	{
+		options.back().sec_message.push_back(msg);
+	}
+	else options[index].sec_message.push_back(msg);
+	return;
+}
+
+void DispHelp::addExample(string ex)
+{
+	examples.push_back(ex);
+	return;
+}
+
+void DispHelp::changeLayerOut(int left,int right)
+{
+	front_space = left; back_space = right;
+	return;
+}
+
+void DispHelp::show()
+{
+	int line_length;
+	string segment,full_message;
+	stringstream ss_message;
+	//获取终端窗口的行列数
+	struct winsize w;
+	ioctl(STDOUT_FILENO, TIOCGWINSZ, &w);
+	//显示头信息
+	full_message = ex_name + " " + version + " - " + descript;
+	ss_message.clear(); ss_message.str(full_message);
+
+	line_length = front_space + back_space;
+	while(ss_message >> segment)
+	{
+		if ((line_length+segment.length()+1) <= w.ws_col)
+		{
+			if (line_length == (front_space + back_space))
+			{
+				for (int i = 0; i < front_space; i++) clog << " ";
+				clog << segment << " ";
+				line_length += (segment.length()+1);
+			}
+			else
+			{
+				clog << segment << " ";
+				line_length += (segment.length()+1);
+			}
+		}
+		else
+		{
+			clog << endl;
+			for (int i = 0; i < front_space+4; i++) clog << " ";
+			clog << segment << " ";
+			line_length = (segment.length()+9+front_space+back_space);
+		}
+	}
+	clog << endl;
+
+	ss_message.clear(); ss_message.str(author);
+	line_length = front_space + back_space;;
+	while(ss_message >> segment)
+	{
+		if ((line_length+segment.length()+1) <= w.ws_col)
+		{
+			if (line_length == (front_space + back_space))
+			{
+				for (int i = 0; i < front_space; i++) clog << " ";
+				clog << "Author: " << segment << " ";
+				line_length += (segment.length()+9);
+			}
+			else
+			{
+				clog << segment << " ";
+				line_length += (segment.length()+1);
+			}
+		}
+		else
+		{
+			clog << endl;
+			for (int i = 0; i < front_space+4; i++) clog << " ";
+			clog << segment << " ";
+			line_length = (segment.length()+9+front_space+back_space);
+		}
+	}
+	clog << endl;
+
+	if (!usages.empty())
+	{
+		for (int i = 0; i < front_space; i++) clog << " ";
+		clog << "Usage:" << endl;
+		for (int i = 0; i < usages.size(); i++)
+		{
+			ss_message.clear(); ss_message.str(usages[i]);
+
+			line_length = front_space + back_space + 4;
+			while(ss_message >> segment)
+			{
+				if ((line_length+segment.length()+1) <= w.ws_col)
+				{
+					if (line_length == (front_space+back_space+4))
+					{
+						for (int i = 0; i < front_space+4; i++) clog << " ";
+						clog << segment << " ";
+						line_length += (segment.length()+1);
+					}
+					else
+					{
+						clog << segment << " ";
+						line_length += (segment.length()+1);
+					}
+				}
+				else
+				{
+					clog << endl;
+					for (int i = 0; i < front_space+9; i++) clog << " ";
+					clog << segment << " ";
+					line_length = (segment.length()+10+front_space+back_space);
+				}
+			}
+			clog << endl;
+		}
+	}
+
+	if (!options.empty())
+	{
+		for (int i = 0; i < front_space; i++) clog << " ";
+		clog << "Options:" << endl;
+		for (int i = 0; i < options.size(); i++)
+		{
+			if (options[i].flag_l == "")
+			{
+				full_message = options[i].flag_s+" : "+options[i].message;
+				ss_message.clear(); ss_message.str(full_message);
+
+				line_length = front_space + back_space + 4;
+				while(ss_message >> segment)
+				{
+					if ((line_length+segment.length()+1) <= w.ws_col)
+					{
+						if (line_length == (front_space+back_space+4))
+						{
+							for (int i = 0; i < front_space+4; i++) clog << " ";
+							clog << segment << " ";
+							line_length += (segment.length()+1);
+						}
+						else
+						{
+							clog << segment << " ";
+							line_length += (segment.length()+1);
+						}
+					}
+					else
+					{
+						clog << endl;
+						for (int i = 0; i < front_space+9; i++) clog << " ";
+						clog << segment << " ";
+						line_length = (segment.length()+10+front_space+back_space);
+					}
+				}
+				clog << endl;
+
+				if (!options[i].sec_message.empty())
+				{
+					for (int j = 0; j < options[i].sec_message.size(); j++)
+					{
+						ss_message.clear(); ss_message.str(options[i].sec_message[j]);
+
+						line_length = front_space + back_space + 9;
+						while(ss_message >> segment)
+						{
+							if ((line_length+segment.length()+1) <= w.ws_col)
+							{
+								if (line_length == (front_space+back_space+9))
+								{
+									for (int i = 0; i < front_space+9; i++) clog << " ";
+									clog << segment << " ";
+									line_length += (segment.length()+1);
+								}
+								else
+								{
+									clog << segment << " ";
+									line_length += (segment.length()+1);
+								}
+							}
+							else
+							{
+								clog << endl;
+								for (int i = 0; i < front_space+13; i++) clog << " ";
+								clog << segment << " ";
+								line_length = (segment.length()+14+front_space+back_space);
+							}
+						}
+						clog << endl;
+					}
+				}
+			}
+			else
+			{
+				full_message = options[i].flag_s+" | "+options[i].flag_l+" : "+options[i].message;
+				ss_message.clear(); ss_message.str(full_message);
+
+				line_length = front_space + back_space + 4;
+				while(ss_message >> segment)
+				{
+					if ((line_length+segment.length()+1) <= w.ws_col)
+					{
+						if (line_length == (front_space+back_space+4))
+						{
+							for (int i = 0; i < front_space+4; i++) clog << " ";
+							clog << segment << " ";
+							line_length += (segment.length()+1);
+						}
+						else
+						{
+							clog << segment << " ";
+							line_length += (segment.length()+1);
+						}
+					}
+					else
+					{
+						clog << endl;
+						for (int i = 0; i < front_space+9; i++) clog << " ";
+						clog << segment << " ";
+						line_length = (segment.length()+10+front_space+back_space);
+					}
+				}
+				clog << endl;
+
+				if (!options[i].sec_message.empty())
+				{
+					for (int j = 0; j < options[i].sec_message.size(); j++)
+					{
+						ss_message.clear(); ss_message.str(options[i].sec_message[j]);
+
+						line_length = front_space + back_space + 9;
+						while(ss_message >> segment)
+						{
+							if ((line_length+segment.length()+1) <= w.ws_col)
+							{
+								if (line_length == (front_space+back_space+9))
+								{
+									for (int i = 0; i < front_space+9; i++) clog << " ";
+									clog << segment << " ";
+									line_length += (segment.length()+1);
+								}
+								else
+								{
+									clog << segment << " ";
+									line_length += (segment.length()+1);
+								}
+							}
+							else
+							{
+								clog << endl;
+								for (int i = 0; i < front_space+13; i++) clog << " ";
+								clog << segment << " ";
+								line_length = (segment.length()+14+front_space+back_space);
+							}
+						}
+						clog << endl;
+					}
+				}
+			}
+		}
+	}
+
+	if (!examples.empty())
+	{
+		for (int i = 0; i < front_space; i++) clog << " ";
+		clog << "Examples:" << endl;
+		for (int i = 0; i < examples.size(); i++)
+		{
+			ss_message.clear(); ss_message.str(examples[i]);
+
+			line_length = front_space + back_space + 4;
+			while(ss_message >> segment)
+			{
+				if ((line_length+segment.length()+1) <= w.ws_col)
+				{
+					if (line_length == (front_space+back_space+4))
+					{
+						for (int i = 0; i < front_space+4; i++) clog << " ";
+						clog << segment << " ";
+						line_length += (segment.length()+1);
+					}
+					else
+					{
+						clog << segment << " ";
+						line_length += (segment.length()+1);
+					}
+				}
+				else
+				{
+					clog << endl;
+					for (int i = 0; i < front_space+9; i++) clog << " ";
+					clog << segment << " ";
+					line_length = (segment.length()+10+front_space+back_space);
+				}
+			}
+			clog << endl;
+		}
+	}
+	return;
+}

src/disp_help.h

@@ -0,0 +1,56 @@
+#ifndef _DISPHELP_H
+#define _DISPHELP_H
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <string.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <unistd.h>
+#include <iomanip>
+#include <sys/ioctl.h>
+#include "vector"
+
+using namespace std;
+
+typedef vector<string> strArray;
+
+struct option
+{
+	string flag_s,flag_l;
+	string message;
+	strArray sec_message;
+	option()
+	{
+		flag_s = flag_l = message = "";
+	}
+};
+typedef vector<option> opArray;
+
+class DispHelp
+{
+public:
+	DispHelp(){
+		front_space = 0;
+		back_space = 10;
+		ex_name = "Execuable";
+		version = "0.0.1";
+		descript = "Brief information about this command.";
+		author = "Author's information.";
+	}
+	~DispHelp(){}
+	void addHeadInfo(string,string,string,string);
+	void addUsage(string);
+	void addOption(string,string,string lflag = "");
+	void addOptionSec(string,int index = -1);
+	void addExample(string);
+	void changeLayerOut(int,int);
+	void show();
+private:
+	string ex_name,version,descript,author;
+	int front_space,back_space;
+	opArray options;
+	strArray examples;
+	strArray usages;
+};
+#endif

src/forward_delta_t.cpp

@@ -0,0 +1,131 @@
+#include "gm3d.h"
+
+int GM3D::ForwardDeltaT(char* noise_level,char* mag_para){
+	int i,j;
+	double I0,A0,I,A;
+	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 noise_mean, noise_dev, temp_noise;
+	if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
+		noise_mean = noise_dev = 0.0;
+	}
+
+	//确定磁化参数
+	if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
+		I0 = I = 90; A0 = A = 0;
+	}
+	else{
+		I0=I0*Pi/180;
+		A0=A0*Pi/180;
+		I=I*Pi/180;
+		A=A*Pi/180;
+
+		k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
+		k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
+		k3=cos(I0)*cos(A0)*cos(I)*sin(A)+cos(I0)*sin(A0)*cos(I)*cos(A);
+		k4=cos(I0)*cos(A0)*cos(I)*cos(A);
+		k5=cos(I0)*sin(A0)*cos(I)*sin(A);
+		k6=-sin(I0)*sin(I);
+	}
+
+	//添加高斯噪声值
+	default_random_engine generator;
+	normal_distribution<double> dist(noise_mean, noise_dev);
+
+	_1dArray temp_obs(model_num_);
+
+	ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
+	for (i = 0; i < obs_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 < model_num_; j++){
+			if (fabs(forward_model_[j]) > ZERO){
+				x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
+				y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
+				z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
+
+				R222=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R122=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R212=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R112=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R221=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R121=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				w222=((x2-obs_p_[i].x)*(x2-obs_p_[i].x))+R222*(z2-obs_p_[i].z)+((z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				w122=((x1-obs_p_[i].x)*(x1-obs_p_[i].x))+R122*(z2-obs_p_[i].z)+((z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				w212=((x2-obs_p_[i].x)*(x2-obs_p_[i].x))+R212*(z2-obs_p_[i].z)+((z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				w112=((x1-obs_p_[i].x)*(x1-obs_p_[i].x))+R112*(z2-obs_p_[i].z)+((z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				w221=((x2-obs_p_[i].x)*(x2-obs_p_[i].x))+R221*(z1-obs_p_[i].z)+((z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				w121=((x1-obs_p_[i].x)*(x1-obs_p_[i].x))+R121*(z1-obs_p_[i].z)+((z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				w211=((x2-obs_p_[i].x)*(x2-obs_p_[i].x))+R211*(z1-obs_p_[i].z)+((z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				w111=((x1-obs_p_[i].x)*(x1-obs_p_[i].x))+R111*(z1-obs_p_[i].z)+((z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				v222=((y2-obs_p_[i].y)*(y2-obs_p_[i].y))+R222*(z2-obs_p_[i].z)+((z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				v122=((y2-obs_p_[i].y)*(y2-obs_p_[i].y))+R122*(z2-obs_p_[i].z)+((z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				v212=((y1-obs_p_[i].y)*(y1-obs_p_[i].y))+R212*(z2-obs_p_[i].z)+((z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				v112=((y1-obs_p_[i].y)*(y1-obs_p_[i].y))+R112*(z2-obs_p_[i].z)+((z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				v221=((y2-obs_p_[i].y)*(y2-obs_p_[i].y))+R221*(z1-obs_p_[i].z)+((z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				v121=((y2-obs_p_[i].y)*(y2-obs_p_[i].y))+R121*(z1-obs_p_[i].z)+((z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				v211=((y1-obs_p_[i].y)*(y1-obs_p_[i].y))+R211*(z1-obs_p_[i].z)+((z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				v111=((y1-obs_p_[i].y)*(y1-obs_p_[i].y))+R111*(z1-obs_p_[i].z)+((z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				G222=k1*log(R222+x2-obs_p_[i].x)+k2*log(R222+y2-obs_p_[i].y)+k3*log(R222+z2-obs_p_[i].z)
+				+k4*atan(((x2-obs_p_[i].x)*(y2-obs_p_[i].y))/w222)+k5*atan(((x2-obs_p_[i].x)*(y2-obs_p_[i].y))/v222)
+				+k6*atan(((x2-obs_p_[i].x)*(y2-obs_p_[i].y))/(R222*(z2-obs_p_[i].z)));
+
+				G122=k1*log(R122+x1-obs_p_[i].x)+k2*log(R122+y2-obs_p_[i].y)+k3*log(R122+z2-obs_p_[i].z)
+				+k4*atan(((x1-obs_p_[i].x)*(y2-obs_p_[i].y))/w122)+k5*atan(((x1-obs_p_[i].x)*(y2-obs_p_[i].y))/v122)
+				+k6*atan(((x1-obs_p_[i].x)*(y2-obs_p_[i].y))/(R122*(z2-obs_p_[i].z)));
+
+				G212=k1*log(R212+x2-obs_p_[i].x)+k2*log(R212+y1-obs_p_[i].y)+k3*log(R212+z2-obs_p_[i].z)
+				+k4*atan(((x2-obs_p_[i].x)*(y1-obs_p_[i].y))/w212)+k5*atan(((x2-obs_p_[i].x)*(y1-obs_p_[i].y))/v212)
+				+k6*atan(((x2-obs_p_[i].x)*(y1-obs_p_[i].y))/(R212*(z2-obs_p_[i].z)));
+
+				G112=k1*log(R112+x1-obs_p_[i].x)+k2*log(R112+y1-obs_p_[i].y)+k3*log(R112+z2-obs_p_[i].z)
+				+k4*atan(((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/w112)+k5*atan(((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/v112)
+				+k6*atan(((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/(R112*(z2-obs_p_[i].z)));
+
+				G221=k1*log(R221+x2-obs_p_[i].x)+k2*log(R221+y2-obs_p_[i].y)+k3*log(R221+z1-obs_p_[i].z)
+				+k4*atan(((x2-obs_p_[i].x)*(y2-obs_p_[i].y))/w221)+k5*atan(((x2-obs_p_[i].x)*(y2-obs_p_[i].y))/v221)
+				+k6*atan(((x2-obs_p_[i].x)*(y2-obs_p_[i].y))/(R221*(z1-obs_p_[i].z)));
+
+				G121=k1*log(R121+x1-obs_p_[i].x)+k2*log(R121+y2-obs_p_[i].y)+k3*log(R121+z1-obs_p_[i].z)
+				+k4*atan(((x1-obs_p_[i].x)*(y2-obs_p_[i].y))/w121)+k5*atan(((x1-obs_p_[i].x)*(y2-obs_p_[i].y))/v121)
+				+k6*atan(((x1-obs_p_[i].x)*(y2-obs_p_[i].y))/(R121*(z1-obs_p_[i].z)));
+
+				G211=k1*log(R211+x2-obs_p_[i].x)+k2*log(R211+y1-obs_p_[i].y)+k3*log(R211+z1-obs_p_[i].z)
+				+k4*atan(((x2-obs_p_[i].x)*(y1-obs_p_[i].y))/w211)+k5*atan(((x2-obs_p_[i].x)*(y1-obs_p_[i].y))/v211)
+				+k6*atan(((x2-obs_p_[i].x)*(y1-obs_p_[i].y))/(R211*(z1-obs_p_[i].z)));
+
+				G111=k1*log(R111+x1-obs_p_[i].x)+k2*log(R111+y1-obs_p_[i].y)+k3*log(R111+z1-obs_p_[i].z)
+				+k4*atan(((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/w111)+k5*atan(((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/v111)
+				+k6*atan(((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/(R111*(z1-obs_p_[i].z)));
+
+				temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
+			}
+		}
+
+		for (int n = 0; n < model_num_; n++){
+			obs_p_[i].val += temp_obs[n];
+		}
+	}
+
+	for (int i = 0; i < obs_num_; i++){
+		temp_noise = dist(generator);
+		obs_p_[i].val += temp_noise;
+		obs_p_[i].dev = noise_dev;
+	}
+	return 0;
+}

src/forward_delta_tx.cpp

@@ -0,0 +1,108 @@
+#include "gm3d.h"
+
+int GM3D::ForwardDeltaTx(char* noise_level,char* mag_para){
+	int i,j;
+	double I0,A0,I,A;
+	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 noise_mean, noise_dev, temp_noise;
+	if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
+		noise_mean = noise_dev = 0.0;
+	}
+
+	//确定磁化参数
+	if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
+		I0 = I = 90; A0 = A = 0;
+	}
+	else{
+		I0=I0*Pi/180;
+		A0=A0*Pi/180;
+		I=I*Pi/180;
+		A=A*Pi/180;
+
+		k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
+		k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
+		k3=cos(I0)*cos(A0)*cos(I)*sin(A)+cos(I0)*sin(A0)*cos(I)*cos(A);
+		k4=cos(I0)*cos(A0)*cos(I)*cos(A);
+		k5=cos(I0)*sin(A0)*cos(I)*sin(A);
+		k6=-sin(I0)*sin(I);
+	}
+
+	//添加高斯噪声值
+	default_random_engine generator;
+	normal_distribution<double> dist(noise_mean, noise_dev);
+
+	_1dArray temp_obs(model_num_);
+
+	ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
+	for (i = 0; i < obs_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 < model_num_; j++){
+			if (fabs(forward_model_[j]) > ZERO){
+				x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
+				y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
+				z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
+
+				R222=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R122=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R212=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R112=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R221=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R121=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				w222=((y2-obs_p_[i].y)*(pow(x2-obs_p_[i].x,2)-R222*(z2-obs_p_[i].z)))/(R222*(pow(x2-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2))*(R222+z2-obs_p_[i].z));
+				w122=((y2-obs_p_[i].y)*(pow(x1-obs_p_[i].x,2)-R122*(z2-obs_p_[i].z)))/(R122*(pow(x1-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2))*(R122+z2-obs_p_[i].z));
+				w212=((y1-obs_p_[i].y)*(pow(x2-obs_p_[i].x,2)-R212*(z2-obs_p_[i].z)))/(R212*(pow(x2-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2))*(R212+z2-obs_p_[i].z));
+				w112=((y1-obs_p_[i].y)*(pow(x1-obs_p_[i].x,2)-R112*(z2-obs_p_[i].z)))/(R112*(pow(x1-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2))*(R112+z2-obs_p_[i].z));
+				w221=((y2-obs_p_[i].y)*(pow(x2-obs_p_[i].x,2)-R221*(z1-obs_p_[i].z)))/(R221*(pow(x2-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2))*(R221+z1-obs_p_[i].z));
+				w121=((y2-obs_p_[i].y)*(pow(x1-obs_p_[i].x,2)-R121*(z1-obs_p_[i].z)))/(R121*(pow(x1-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2))*(R121+z1-obs_p_[i].z));
+				w211=((y1-obs_p_[i].y)*(pow(x2-obs_p_[i].x,2)-R211*(z1-obs_p_[i].z)))/(R211*(pow(x2-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2))*(R211+z1-obs_p_[i].z));
+				w111=((y1-obs_p_[i].y)*(pow(x1-obs_p_[i].x,2)-R111*(z1-obs_p_[i].z)))/(R111*(pow(x1-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2))*(R111+z1-obs_p_[i].z));
+
+				v222=((y2-obs_p_[i].y)*(z2-obs_p_[i].z))/(R222*(pow(x2-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				v122=((y2-obs_p_[i].y)*(z2-obs_p_[i].z))/(R122*(pow(x1-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				v212=((y1-obs_p_[i].y)*(z2-obs_p_[i].z))/(R212*(pow(x2-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				v112=((y1-obs_p_[i].y)*(z2-obs_p_[i].z))/(R112*(pow(x1-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				v221=((y2-obs_p_[i].y)*(z1-obs_p_[i].z))/(R221*(pow(x2-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+				v121=((y2-obs_p_[i].y)*(z1-obs_p_[i].z))/(R121*(pow(x1-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+				v211=((y1-obs_p_[i].y)*(z1-obs_p_[i].z))/(R211*(pow(x2-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+				v111=((y1-obs_p_[i].y)*(z1-obs_p_[i].z))/(R111*(pow(x1-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+
+				G222=-k1*pow(R222,-1)-k2*(x2-obs_p_[i].x)/(R222*(R222+y2-obs_p_[i].y))-k3*(x2-obs_p_[i].x)/(R222*(R222+z2-obs_p_[i].z))+k4*w222-k5*(y2-obs_p_[i].y)/(R222*(R222+z2-obs_p_[i].z))-k6*v222;
+				G122=-k1*pow(R122,-1)-k2*(x1-obs_p_[i].x)/(R122*(R122+y2-obs_p_[i].y))-k3*(x1-obs_p_[i].x)/(R122*(R122+z2-obs_p_[i].z))+k4*w122-k5*(y2-obs_p_[i].y)/(R122*(R122+z2-obs_p_[i].z))-k6*v122;
+				G212=-k1*pow(R212,-1)-k2*(x2-obs_p_[i].x)/(R212*(R212+y1-obs_p_[i].y))-k3*(x2-obs_p_[i].x)/(R212*(R212+z2-obs_p_[i].z))+k4*w212-k5*(y1-obs_p_[i].y)/(R212*(R212+z2-obs_p_[i].z))-k6*v212;
+				G112=-k1*pow(R112,-1)-k2*(x1-obs_p_[i].x)/(R112*(R112+y1-obs_p_[i].y))-k3*(x1-obs_p_[i].x)/(R112*(R112+z2-obs_p_[i].z))+k4*w112-k5*(y1-obs_p_[i].y)/(R112*(R112+z2-obs_p_[i].z))-k6*v112;
+				G221=-k1*pow(R221,-1)-k2*(x2-obs_p_[i].x)/(R221*(R221+y2-obs_p_[i].y))-k3*(x2-obs_p_[i].x)/(R221*(R221+z1-obs_p_[i].z))+k4*w221-k5*(y2-obs_p_[i].y)/(R221*(R221+z1-obs_p_[i].z))-k6*v221;
+				G121=-k1*pow(R121,-1)-k2*(x1-obs_p_[i].x)/(R121*(R121+y2-obs_p_[i].y))-k3*(x1-obs_p_[i].x)/(R121*(R121+z1-obs_p_[i].z))+k4*w121-k5*(y2-obs_p_[i].y)/(R121*(R121+z1-obs_p_[i].z))-k6*v121;
+				G211=-k1*pow(R211,-1)-k2*(x2-obs_p_[i].x)/(R211*(R211+y1-obs_p_[i].y))-k3*(x2-obs_p_[i].x)/(R211*(R211+z1-obs_p_[i].z))+k4*w211-k5*(y1-obs_p_[i].y)/(R211*(R211+z1-obs_p_[i].z))-k6*v211;
+				G111=-k1*pow(R111,-1)-k2*(x1-obs_p_[i].x)/(R111*(R111+y1-obs_p_[i].y))-k3*(x1-obs_p_[i].x)/(R111*(R111+z1-obs_p_[i].z))+k4*w111-k5*(y1-obs_p_[i].y)/(R111*(R111+z1-obs_p_[i].z))-k6*v111;
+
+				temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
+			}
+		}
+
+		for (int n = 0; n < model_num_; n++){
+			obs_p_[i].val += temp_obs[n];
+		}
+	}
+
+	for (int i = 0; i < obs_num_; i++){
+		temp_noise = dist(generator);
+		obs_p_[i].val += temp_noise;
+		obs_p_[i].dev = noise_dev;
+	}
+	return 0;
+}

src/forward_delta_ty.cpp

@@ -0,0 +1,108 @@
+#include "gm3d.h"
+
+int GM3D::ForwardDeltaTy(char* noise_level,char* mag_para){
+	int i,j;
+	double I0,A0,I,A;
+	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 noise_mean, noise_dev, temp_noise;
+	if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
+		noise_mean = noise_dev = 0.0;
+	}
+
+	//确定磁化参数
+	if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
+		I0 = I = 90; A0 = A = 0;
+	}
+	else{
+		I0=I0*Pi/180;
+		A0=A0*Pi/180;
+		I=I*Pi/180;
+		A=A*Pi/180;
+
+		k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
+		k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
+		k3=cos(I0)*cos(A0)*cos(I)*sin(A)+cos(I0)*sin(A0)*cos(I)*cos(A);
+		k4=cos(I0)*cos(A0)*cos(I)*cos(A);
+		k5=cos(I0)*sin(A0)*cos(I)*sin(A);
+		k6=-sin(I0)*sin(I);
+	}
+
+	//添加高斯噪声值
+	default_random_engine generator;
+	normal_distribution<double> dist(noise_mean, noise_dev);
+
+	_1dArray temp_obs(model_num_);
+
+	ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
+	for (i = 0; i < obs_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 < model_num_; j++){
+			if (fabs(forward_model_[j]) > ZERO){
+				x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
+				y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
+				z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
+
+				R222=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R122=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R212=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R112=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R221=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R121=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				w222=((x2-obs_p_[i].x)*(pow(y2-obs_p_[i].y,2)-R222*(z2-obs_p_[i].z)))/(R222*(pow(y2-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))*(R222+z2-obs_p_[i].z));
+				w122=((x1-obs_p_[i].x)*(pow(y2-obs_p_[i].y,2)-R122*(z2-obs_p_[i].z)))/(R122*(pow(y2-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))*(R122+z2-obs_p_[i].z));
+				w212=((x2-obs_p_[i].x)*(pow(y1-obs_p_[i].y,2)-R212*(z2-obs_p_[i].z)))/(R212*(pow(y1-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))*(R212+z2-obs_p_[i].z));
+				w112=((x1-obs_p_[i].x)*(pow(y1-obs_p_[i].y,2)-R112*(z2-obs_p_[i].z)))/(R112*(pow(y1-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))*(R112+z2-obs_p_[i].z));
+				w221=((x2-obs_p_[i].x)*(pow(y2-obs_p_[i].y,2)-R221*(z1-obs_p_[i].z)))/(R221*(pow(y2-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))*(R221+z1-obs_p_[i].z));
+				w121=((x1-obs_p_[i].x)*(pow(y2-obs_p_[i].y,2)-R121*(z1-obs_p_[i].z)))/(R121*(pow(y2-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))*(R121+z1-obs_p_[i].z));
+				w211=((x2-obs_p_[i].x)*(pow(y1-obs_p_[i].y,2)-R211*(z1-obs_p_[i].z)))/(R211*(pow(y1-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))*(R211+z1-obs_p_[i].z));
+				w111=((x1-obs_p_[i].x)*(pow(y1-obs_p_[i].y,2)-R111*(z1-obs_p_[i].z)))/(R111*(pow(y1-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))*(R111+z1-obs_p_[i].z));
+
+				v222=((x2-obs_p_[i].x)*(z2-obs_p_[i].z))/(R222*(pow(y2-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2)));
+				v122=((x1-obs_p_[i].x)*(z2-obs_p_[i].z))/(R122*(pow(y2-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2)));
+				v212=((x2-obs_p_[i].x)*(z2-obs_p_[i].z))/(R212*(pow(y1-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2)));
+				v112=((x1-obs_p_[i].x)*(z2-obs_p_[i].z))/(R112*(pow(y1-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2)));
+				v221=((x2-obs_p_[i].x)*(z1-obs_p_[i].z))/(R221*(pow(y2-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2)));
+				v121=((x1-obs_p_[i].x)*(z1-obs_p_[i].z))/(R121*(pow(y2-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2)));
+				v211=((x2-obs_p_[i].x)*(z1-obs_p_[i].z))/(R211*(pow(y1-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2)));
+				v111=((x1-obs_p_[i].x)*(z1-obs_p_[i].z))/(R111*(pow(y1-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2)));
+				
+				G222=-k1*(y2-obs_p_[i].y)/(R222*(R222+x2-obs_p_[i].x))-k2*pow(R222,-1)-k3*(y2-obs_p_[i].y)/(R222*(R222+z2-obs_p_[i].z))-k4*(x2-obs_p_[i].x)/(R222*(R222+z2-obs_p_[i].z))+k5*w222-k6*v222;
+				G122=-k1*(y2-obs_p_[i].y)/(R122*(R122+x1-obs_p_[i].x))-k2*pow(R122,-1)-k3*(y2-obs_p_[i].y)/(R122*(R122+z2-obs_p_[i].z))-k4*(x1-obs_p_[i].x)/(R122*(R122+z2-obs_p_[i].z))+k5*w122-k6*v122;
+				G212=-k1*(y1-obs_p_[i].y)/(R212*(R212+x2-obs_p_[i].x))-k2*pow(R212,-1)-k3*(y1-obs_p_[i].y)/(R212*(R212+z2-obs_p_[i].z))-k4*(x2-obs_p_[i].x)/(R212*(R212+z2-obs_p_[i].z))+k5*w212-k6*v212;
+				G112=-k1*(y1-obs_p_[i].y)/(R112*(R112+x1-obs_p_[i].x))-k2*pow(R112,-1)-k3*(y1-obs_p_[i].y)/(R112*(R112+z2-obs_p_[i].z))-k4*(x1-obs_p_[i].x)/(R112*(R112+z2-obs_p_[i].z))+k5*w112-k6*v112;
+				G221=-k1*(y2-obs_p_[i].y)/(R221*(R221+x2-obs_p_[i].x))-k2*pow(R221,-1)-k3*(y2-obs_p_[i].y)/(R221*(R221+z1-obs_p_[i].z))-k4*(x2-obs_p_[i].x)/(R221*(R221+z1-obs_p_[i].z))+k5*w221-k6*v221;
+				G121=-k1*(y2-obs_p_[i].y)/(R121*(R121+x1-obs_p_[i].x))-k2*pow(R121,-1)-k3*(y2-obs_p_[i].y)/(R121*(R121+z1-obs_p_[i].z))-k4*(x1-obs_p_[i].x)/(R121*(R121+z1-obs_p_[i].z))+k5*w121-k6*v121;
+				G211=-k1*(y1-obs_p_[i].y)/(R211*(R211+x2-obs_p_[i].x))-k2*pow(R211,-1)-k3*(y1-obs_p_[i].y)/(R211*(R211+z1-obs_p_[i].z))-k4*(x2-obs_p_[i].x)/(R211*(R211+z1-obs_p_[i].z))+k5*w211-k6*v211;
+				G111=-k1*(y1-obs_p_[i].y)/(R111*(R111+x1-obs_p_[i].x))-k2*pow(R111,-1)-k3*(y1-obs_p_[i].y)/(R111*(R111+z1-obs_p_[i].z))-k4*(x1-obs_p_[i].x)/(R111*(R111+z1-obs_p_[i].z))+k5*w111-k6*v111;
+
+				temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
+			}
+		}
+
+		for (int n = 0; n < model_num_; n++){
+			obs_p_[i].val += temp_obs[n];
+		}
+	}
+
+	for (int i = 0; i < obs_num_; i++){
+		temp_noise = dist(generator);
+		obs_p_[i].val += temp_noise;
+		obs_p_[i].dev = noise_dev;
+	}
+	return 0;
+}

src/forward_delta_tz.cpp

@@ -0,0 +1,119 @@
+#include "gm3d.h"
+
+int GM3D::ForwardDeltaTz(char* noise_level,char* mag_para){
+	int i,j;
+	double I0,A0,I,A;
+	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 t111,t112,t121,t122,t211,t212,t221,t222;
+	double R222,R122,R212,R112,R221,R121,R211,R111;
+	double G222,G122,G212,G112,G221,G121,G211,G111;
+
+	//确定噪声水平
+	double noise_mean, noise_dev, temp_noise;
+	if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
+		noise_mean = noise_dev = 0.0;
+	}
+
+	//确定磁化参数
+	if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
+		I0 = I = 90; A0 = A = 0;
+	}
+	else{
+		I0=I0*Pi/180;
+		A0=A0*Pi/180;
+		I=I*Pi/180;
+		A=A*Pi/180;
+
+		k1=cos(I0)*sin(A0)*sin(I)+sin(I0)*cos(I)*sin(A);
+		k2=cos(I0)*cos(A0)*sin(I)+sin(I0)*cos(I)*cos(A);
+		k3=cos(I0)*cos(A0)*cos(I)*sin(A)+cos(I0)*sin(A0)*cos(I)*cos(A);
+		k4=cos(I0)*cos(A0)*cos(I)*cos(A);
+		k5=cos(I0)*sin(A0)*cos(I)*sin(A);
+		k6=-sin(I0)*sin(I);
+	}
+
+	//添加高斯噪声值
+	default_random_engine generator;
+	normal_distribution<double> dist(noise_mean, noise_dev);
+
+	_1dArray temp_obs(model_num_);
+
+	ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
+	for (i = 0; i < obs_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,\
+		t111,t112,t121,t122,t211,t212,t221,t222,\
+		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 < model_num_; j++){
+			if (fabs(forward_model_[j]) > ZERO){
+				x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
+				y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
+				z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
+
+				R222=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R122=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R212=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R112=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R221=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R121=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				w222=((x2-obs_p_[i].x)*(y2-obs_p_[i].y)*(pow(z2-obs_p_[i].z,2)+pow(R222,2)))/(R222*(pow(y2-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))*(pow(x2-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				w122=((x1-obs_p_[i].x)*(y2-obs_p_[i].y)*(pow(z2-obs_p_[i].z,2)+pow(R122,2)))/(R122*(pow(y2-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))*(pow(x1-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				w212=((x2-obs_p_[i].x)*(y1-obs_p_[i].y)*(pow(z2-obs_p_[i].z,2)+pow(R212,2)))/(R212*(pow(y1-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))*(pow(x2-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				w112=((x1-obs_p_[i].x)*(y1-obs_p_[i].y)*(pow(z2-obs_p_[i].z,2)+pow(R112,2)))/(R112*(pow(y1-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))*(pow(x1-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				w221=((x2-obs_p_[i].x)*(y2-obs_p_[i].y)*(pow(z1-obs_p_[i].z,2)+pow(R221,2)))/(R221*(pow(y2-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))*(pow(x2-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+				w121=((x1-obs_p_[i].x)*(y2-obs_p_[i].y)*(pow(z1-obs_p_[i].z,2)+pow(R121,2)))/(R121*(pow(y2-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))*(pow(x1-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+				w211=((x2-obs_p_[i].x)*(y1-obs_p_[i].y)*(pow(z1-obs_p_[i].z,2)+pow(R211,2)))/(R211*(pow(y1-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))*(pow(x2-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+				w111=((x1-obs_p_[i].x)*(y1-obs_p_[i].y)*(pow(z1-obs_p_[i].z,2)+pow(R111,2)))/(R111*(pow(y1-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))*(pow(x1-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+
+				v222=((x2-obs_p_[i].x)*(y2-obs_p_[i].y))/(R222*(pow(y2-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2)));
+				v122=((x1-obs_p_[i].x)*(y2-obs_p_[i].y))/(R122*(pow(y2-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2)));
+				v212=((x2-obs_p_[i].x)*(y1-obs_p_[i].y))/(R212*(pow(y1-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2)));
+				v112=((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/(R112*(pow(y1-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2)));
+				v221=((x2-obs_p_[i].x)*(y2-obs_p_[i].y))/(R221*(pow(y2-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2)));
+				v121=((x1-obs_p_[i].x)*(y2-obs_p_[i].y))/(R121*(pow(y2-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2)));
+				v211=((x2-obs_p_[i].x)*(y1-obs_p_[i].y))/(R211*(pow(y1-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2)));
+				v111=((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/(R111*(pow(y1-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2)));
+
+				t222=((x2-obs_p_[i].x)*(y2-obs_p_[i].y))/(R222*(pow(x2-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				t122=((x1-obs_p_[i].x)*(y2-obs_p_[i].y))/(R122*(pow(x1-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				t212=((x2-obs_p_[i].x)*(y1-obs_p_[i].y))/(R212*(pow(x2-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				t112=((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/(R112*(pow(x1-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2)));
+				t221=((x2-obs_p_[i].x)*(y2-obs_p_[i].y))/(R221*(pow(x2-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+				t121=((x1-obs_p_[i].x)*(y2-obs_p_[i].y))/(R121*(pow(x1-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+				t211=((x2-obs_p_[i].x)*(y1-obs_p_[i].y))/(R211*(pow(x2-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+				t111=((x1-obs_p_[i].x)*(y1-obs_p_[i].y))/(R111*(pow(x1-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2)));
+				
+				G222=-k1*(z2-obs_p_[i].z)/(R222*(R222+x2-obs_p_[i].x))-k2*(z2-obs_p_[i].z)/(R222*(R222+y2-obs_p_[i].y))-k3/R222+k4*t222+k5*v222+k6*w222;
+				G122=-k1*(z2-obs_p_[i].z)/(R122*(R122+x1-obs_p_[i].x))-k2*(z2-obs_p_[i].z)/(R122*(R122+y2-obs_p_[i].y))-k3/R122+k4*t122+k5*v122+k6*w122;
+				G212=-k1*(z2-obs_p_[i].z)/(R212*(R212+x2-obs_p_[i].x))-k2*(z2-obs_p_[i].z)/(R212*(R212+y1-obs_p_[i].y))-k3/R212+k4*t212+k5*v212+k6*w212;
+				G112=-k1*(z2-obs_p_[i].z)/(R112*(R112+x1-obs_p_[i].x))-k2*(z2-obs_p_[i].z)/(R112*(R112+y1-obs_p_[i].y))-k3/R112+k4*t112+k5*v112+k6*w112;
+				G221=-k1*(z1-obs_p_[i].z)/(R221*(R221+x2-obs_p_[i].x))-k2*(z1-obs_p_[i].z)/(R221*(R221+y2-obs_p_[i].y))-k3/R221+k4*t221+k5*v221+k6*w221;
+				G121=-k1*(z1-obs_p_[i].z)/(R121*(R121+x1-obs_p_[i].x))-k2*(z1-obs_p_[i].z)/(R121*(R121+y2-obs_p_[i].y))-k3/R121+k4*t121+k5*v121+k6*w121;
+				G211=-k1*(z1-obs_p_[i].z)/(R211*(R211+x2-obs_p_[i].x))-k2*(z1-obs_p_[i].z)/(R211*(R211+y1-obs_p_[i].y))-k3/R211+k4*t211+k5*v211+k6*w211;
+				G111=-k1*(z1-obs_p_[i].z)/(R111*(R111+x1-obs_p_[i].x))-k2*(z1-obs_p_[i].z)/(R111*(R111+y1-obs_p_[i].y))-k3/R111+k4*t111+k5*v111+k6*w111;
+
+				temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
+			}
+		}
+
+		for (int n = 0; n < model_num_; n++){
+			obs_p_[i].val += temp_obs[n];
+		}
+	}
+
+	for (int i = 0; i < obs_num_; i++){
+		temp_noise = dist(generator);
+		obs_p_[i].val += temp_noise;
+		obs_p_[i].dev = noise_dev;
+	}
+	return 0;
+}

src/forward_hax.cpp

@@ -0,0 +1,103 @@
+#include "gm3d.h"
+
+int GM3D::ForwardHax(char* noise_level,char* mag_para){
+	int i,j;
+	double I0,A0,I,A;
+	double Alpha,Beta,Gamma;
+	double x1,x2,y1,y2,z1,z2;
+	double w111_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1;
+	double w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2;
+	double R222,R122,R212,R112,R221,R121,R211,R111;
+	double G222,G122,G212,G112,G221,G121,G211,G111;
+
+	//确定噪声水平
+	double noise_mean, noise_dev, temp_noise;
+	if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
+		noise_mean = noise_dev = 0.0;
+	}
+
+	//确定磁化参数
+	if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
+		I = 90; A = 0;
+	}
+	else{
+		I=I*Pi/180;
+		A=A*Pi/180;
+
+		Alpha=cos(I)*cos(A);
+		Beta=cos(I)*sin(A);
+		Gamma=sin(I);
+	}
+
+	//添加高斯噪声值
+	default_random_engine generator;
+	normal_distribution<double> dist(noise_mean, noise_dev);
+
+	_1dArray temp_obs(model_num_);
+
+	ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
+	for (i = 0; i < obs_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_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1,\
+		w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2,\
+		x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
+		for (j = 0; j < model_num_; j++){
+			if (fabs(forward_model_[j]) > ZERO){
+				x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
+				y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
+				z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
+
+				R222=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R122=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R212=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R112=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R221=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R121=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				w222_1=(x2-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w122_1=(x1-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w212_1=(x2-obs_p_[i].x)*(y1-obs_p_[i].y);
+				w112_1=(x1-obs_p_[i].x)*(y1-obs_p_[i].y);
+				w221_1=(x2-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w121_1=(x1-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w211_1=(x2-obs_p_[i].x)*(y1-obs_p_[i].y);
+				w111_1=(x1-obs_p_[i].x)*(y1-obs_p_[i].y);
+
+				w222_2=(pow(x2-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2))+R222*(z2-obs_p_[i].z);
+				w122_2=(pow(x1-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2))+R122*(z2-obs_p_[i].z);
+				w212_2=(pow(x2-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2))+R212*(z2-obs_p_[i].z);
+				w112_2=(pow(x1-obs_p_[i].x,2)+pow(z2-obs_p_[i].z,2))+R112*(z2-obs_p_[i].z);
+				w221_2=(pow(x2-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2))+R221*(z1-obs_p_[i].z);
+				w121_2=(pow(x1-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2))+R121*(z1-obs_p_[i].z);
+				w211_2=(pow(x2-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2))+R211*(z1-obs_p_[i].z);
+				w111_2=(pow(x1-obs_p_[i].x,2)+pow(z1-obs_p_[i].z,2))+R111*(z1-obs_p_[i].z);
+				
+				G222=Alpha*atan2(w222_1,w222_2)+Beta*log(R222+z2-obs_p_[i].z)+Gamma*log(R222+y2-obs_p_[i].y);
+				G122=Alpha*atan2(w122_1,w122_2)+Beta*log(R122+z2-obs_p_[i].z)+Gamma*log(R122+y2-obs_p_[i].y);
+				G212=Alpha*atan2(w212_1,w212_2)+Beta*log(R212+z2-obs_p_[i].z)+Gamma*log(R212+y1-obs_p_[i].y);
+				G112=Alpha*atan2(w112_1,w112_2)+Beta*log(R112+z2-obs_p_[i].z)+Gamma*log(R112+y1-obs_p_[i].y);
+				G221=Alpha*atan2(w221_1,w221_2)+Beta*log(R221+z1-obs_p_[i].z)+Gamma*log(R221+y2-obs_p_[i].y);
+				G121=Alpha*atan2(w121_1,w121_2)+Beta*log(R121+z1-obs_p_[i].z)+Gamma*log(R121+y2-obs_p_[i].y);
+				G211=Alpha*atan2(w211_1,w211_2)+Beta*log(R211+z1-obs_p_[i].z)+Gamma*log(R211+y1-obs_p_[i].y);
+				G111=Alpha*atan2(w111_1,w111_2)+Beta*log(R111+z1-obs_p_[i].z)+Gamma*log(R111+y1-obs_p_[i].y);
+
+				temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
+			}
+		}
+
+		for (int n = 0; n < model_num_; n++){
+			obs_p_[i].val += temp_obs[n];
+		}
+	}
+
+	for (int i = 0; i < obs_num_; i++){
+		temp_noise = dist(generator);
+		obs_p_[i].val += temp_noise;
+		obs_p_[i].dev = noise_dev;
+	}
+	return 0;
+}

src/forward_hay.cpp

@@ -0,0 +1,103 @@
+#include "gm3d.h"
+
+int GM3D::ForwardHay(char* noise_level,char* mag_para){
+	int i,j;
+	double I0,A0,I,A;
+	double Alpha,Beta,Gamma;
+	double x1,x2,y1,y2,z1,z2;
+	double w111_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1;
+	double w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2;
+	double R222,R122,R212,R112,R221,R121,R211,R111;
+	double G222,G122,G212,G112,G221,G121,G211,G111;
+
+	//确定噪声水平
+	double noise_mean, noise_dev, temp_noise;
+	if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
+		noise_mean = noise_dev = 0.0;
+	}
+
+	//确定磁化参数
+	if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
+		I = 90; A = 0;
+	}
+	else{
+		I=I*Pi/180;
+		A=A*Pi/180;
+
+		Alpha=cos(I)*cos(A);
+		Beta=cos(I)*sin(A);
+		Gamma=sin(I);
+	}
+
+	//添加高斯噪声值
+	default_random_engine generator;
+	normal_distribution<double> dist(noise_mean, noise_dev);
+
+	_1dArray temp_obs(model_num_);
+
+	ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
+	for (i = 0; i < obs_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_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1,\
+		w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2,\
+		x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
+		for (j = 0; j < model_num_; j++){
+			if (fabs(forward_model_[j]) > ZERO){
+				x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
+				y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
+				z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
+
+				R222=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R122=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R212=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R112=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R221=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R121=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				w222_1=(x2-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w122_1=(x1-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w212_1=(x2-obs_p_[i].x)*(y1-obs_p_[i].y);
+				w112_1=(x1-obs_p_[i].x)*(y1-obs_p_[i].y);
+				w221_1=(x2-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w121_1=(x1-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w211_1=(x2-obs_p_[i].x)*(y1-obs_p_[i].y);
+				w111_1=(x1-obs_p_[i].x)*(y1-obs_p_[i].y);
+
+				w222_2=(pow(y2-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))+R222*(z2-obs_p_[i].z);
+				w122_2=(pow(y2-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))+R122*(z2-obs_p_[i].z);
+				w212_2=(pow(y1-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))+R212*(z2-obs_p_[i].z);
+				w112_2=(pow(y1-obs_p_[i].y,2)+pow(z2-obs_p_[i].z,2))+R112*(z2-obs_p_[i].z);
+				w221_2=(pow(y2-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))+R221*(z1-obs_p_[i].z);
+				w121_2=(pow(y2-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))+R121*(z1-obs_p_[i].z);
+				w211_2=(pow(y1-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))+R211*(z1-obs_p_[i].z);
+				w111_2=(pow(y1-obs_p_[i].y,2)+pow(z1-obs_p_[i].z,2))+R111*(z1-obs_p_[i].z);
+				
+				G222=Beta*atan2(w222_1,w222_2)+Alpha*log(R222+z2-obs_p_[i].z)+Gamma*log(R222+x2-obs_p_[i].x);
+				G122=Beta*atan2(w122_1,w122_2)+Alpha*log(R122+z2-obs_p_[i].z)+Gamma*log(R122+x1-obs_p_[i].x);
+				G212=Beta*atan2(w212_1,w212_2)+Alpha*log(R212+z2-obs_p_[i].z)+Gamma*log(R212+x2-obs_p_[i].x);
+				G112=Beta*atan2(w112_1,w112_2)+Alpha*log(R112+z2-obs_p_[i].z)+Gamma*log(R112+x1-obs_p_[i].x);
+				G221=Beta*atan2(w221_1,w221_2)+Alpha*log(R221+z1-obs_p_[i].z)+Gamma*log(R221+x2-obs_p_[i].x);
+				G121=Beta*atan2(w121_1,w121_2)+Alpha*log(R121+z1-obs_p_[i].z)+Gamma*log(R121+x1-obs_p_[i].x);
+				G211=Beta*atan2(w211_1,w211_2)+Alpha*log(R211+z1-obs_p_[i].z)+Gamma*log(R211+x2-obs_p_[i].x);
+				G111=Beta*atan2(w111_1,w111_2)+Alpha*log(R111+z1-obs_p_[i].z)+Gamma*log(R111+x1-obs_p_[i].x);
+
+				temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
+			}
+		}
+
+		for (int n = 0; n < model_num_; n++){
+			obs_p_[i].val += temp_obs[n];
+		}
+	}
+
+	for (int i = 0; i < obs_num_; i++){
+		temp_noise = dist(generator);
+		obs_p_[i].val += temp_noise;
+		obs_p_[i].dev = noise_dev;
+	}
+	return 0;
+}

src/forward_vz.cpp

@@ -0,0 +1,64 @@
+#include "gm3d.h"
+
+int GM3D::ForwardVz(char* noise_level){
+	int i,j;
+	double x1,x2,y1,y2,z1,z2;
+	double R222,R122,R212,R112,R221,R121,R211,R111;
+	double G222,G122,G212,G112,G221,G121,G211,G111;
+
+	//确定噪声水平
+	double noise_mean, noise_dev, temp_noise;
+	if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
+		noise_mean = noise_dev = 0.0;
+	}
+
+	//添加高斯噪声值
+	default_random_engine generator;
+	normal_distribution<double> dist(noise_mean, noise_dev);
+
+	_1dArray temp_obs(model_num_);
+
+	ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
+	for (i = 0; i < obs_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,x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
+		for (j = 0; j < model_num_; j++){
+			if (fabs(forward_model_[j]) > ZERO){
+				x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
+				y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
+				z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
+
+				R222=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R122=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R212=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R112=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R221=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R121=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				G222=(x2-obs_p_[i].x)*log((y2-obs_p_[i].y)+R222)+(y2-obs_p_[i].y)*log((x2-obs_p_[i].x)+R222)+(z2-obs_p_[i].z)*arctg((z2-obs_p_[i].z)*R222/(x2-obs_p_[i].x)/(y2-obs_p_[i].y));
+				G122=(x1-obs_p_[i].x)*log((y2-obs_p_[i].y)+R122)+(y2-obs_p_[i].y)*log((x1-obs_p_[i].x)+R122)+(z2-obs_p_[i].z)*arctg((z2-obs_p_[i].z)*R122/(x1-obs_p_[i].x)/(y2-obs_p_[i].y));
+				G212=(x2-obs_p_[i].x)*log((y1-obs_p_[i].y)+R212)+(y1-obs_p_[i].y)*log((x2-obs_p_[i].x)+R212)+(z2-obs_p_[i].z)*arctg((z2-obs_p_[i].z)*R212/(x2-obs_p_[i].x)/(y1-obs_p_[i].y));
+				G112=(x1-obs_p_[i].x)*log((y1-obs_p_[i].y)+R112)+(y1-obs_p_[i].y)*log((x1-obs_p_[i].x)+R112)+(z2-obs_p_[i].z)*arctg((z2-obs_p_[i].z)*R112/(x1-obs_p_[i].x)/(y1-obs_p_[i].y));
+				G221=(x2-obs_p_[i].x)*log((y2-obs_p_[i].y)+R221)+(y2-obs_p_[i].y)*log((x2-obs_p_[i].x)+R221)+(z1-obs_p_[i].z)*arctg((z1-obs_p_[i].z)*R221/(x2-obs_p_[i].x)/(y2-obs_p_[i].y));
+				G121=(x1-obs_p_[i].x)*log((y2-obs_p_[i].y)+R121)+(y2-obs_p_[i].y)*log((x1-obs_p_[i].x)+R121)+(z1-obs_p_[i].z)*arctg((z1-obs_p_[i].z)*R121/(x1-obs_p_[i].x)/(y2-obs_p_[i].y));
+				G211=(x2-obs_p_[i].x)*log((y1-obs_p_[i].y)+R211)+(y1-obs_p_[i].y)*log((x2-obs_p_[i].x)+R211)+(z1-obs_p_[i].z)*arctg((z1-obs_p_[i].z)*R211/(x2-obs_p_[i].x)/(y1-obs_p_[i].y));
+				G111=(x1-obs_p_[i].x)*log((y1-obs_p_[i].y)+R111)+(y1-obs_p_[i].y)*log((x1-obs_p_[i].x)+R111)+(z1-obs_p_[i].z)*arctg((z1-obs_p_[i].z)*R111/(x1-obs_p_[i].x)/(y1-obs_p_[i].y));
+
+				temp_obs[j] = -1.0*G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
+			}
+		}
+
+		for (int n = 0; n < model_num_; n++){
+			obs_p_[i].val += temp_obs[n];
+		}
+	}
+
+	for (int i = 0; i < obs_num_; i++){
+		temp_noise = dist(generator);
+		obs_p_[i].val += temp_noise;
+		obs_p_[i].dev = noise_dev;
+	}
+	return 0;
+}

src/forward_vzx.cpp

@@ -0,0 +1,63 @@
+#include "gm3d.h"
+
+int GM3D::ForwardVzx(char* noise_level){
+	int i,j;
+	double x1,x2,y1,y2,z1,z2;
+	double R222,R122,R212,R112,R221,R121,R211,R111;
+	double G222,G122,G212,G112,G221,G121,G211,G111;
+
+	double noise_mean, noise_dev, temp_noise;
+	if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
+		noise_mean = noise_dev = 0.0;
+	}
+
+	//添加高斯噪声值
+	default_random_engine generator;
+	normal_distribution<double> dist(noise_mean, noise_dev);
+
+	_1dArray temp_obs(model_num_);
+
+	ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
+	for (i = 0; i < obs_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,x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
+		for (j = 0; j < model_num_; j++){
+			if (fabs(forward_model_[j]) > ZERO){
+				x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
+				y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
+				z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
+
+				R222=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R122=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R212=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R112=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R221=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R121=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				G222=log((y2-obs_p_[i].y)+R222);
+				G122=log((y2-obs_p_[i].y)+R122);
+				G212=log((y1-obs_p_[i].y)+R212);
+				G112=log((y1-obs_p_[i].y)+R112);
+				G221=log((y2-obs_p_[i].y)+R221);
+				G121=log((y2-obs_p_[i].y)+R121);
+				G211=log((y1-obs_p_[i].y)+R211);
+				G111=log((y1-obs_p_[i].y)+R111);
+
+				temp_obs[j] = G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
+			}
+		}
+
+		for (int n = 0; n < model_num_; n++){
+			obs_p_[i].val += temp_obs[n];
+		}
+	}
+
+	for (int i = 0; i < obs_num_; i++){
+		temp_noise = dist(generator);
+		obs_p_[i].val += temp_noise;
+		obs_p_[i].dev = noise_dev;
+	}
+	return 0;
+}

src/forward_vzy.cpp

@@ -0,0 +1,63 @@
+#include "gm3d.h"
+
+int GM3D::ForwardVzy(char* noise_level){
+	int i,j;
+	double x1,x2,y1,y2,z1,z2;
+	double R222,R122,R212,R112,R221,R121,R211,R111;
+	double G222,G122,G212,G112,G221,G121,G211,G111;
+
+	double noise_mean, noise_dev, temp_noise;
+	if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
+		noise_mean = noise_dev = 0.0;
+	}
+
+	//添加高斯噪声值
+	default_random_engine generator;
+	normal_distribution<double> dist(noise_mean, noise_dev);
+
+	_1dArray temp_obs(model_num_);
+
+	ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
+	for (i = 0; i < obs_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,x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
+		for (j = 0; j < model_num_; j++){
+			if (fabs(forward_model_[j]) > ZERO){
+				x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
+				y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
+				z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
+
+				R222=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R122=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R212=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R112=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R221=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R121=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				G222=log((x2-obs_p_[i].x)+R222);
+				G122=log((x1-obs_p_[i].x)+R122);
+				G212=log((x2-obs_p_[i].x)+R212);
+				G112=log((x1-obs_p_[i].x)+R112);
+				G221=log((x2-obs_p_[i].x)+R221);
+				G121=log((x1-obs_p_[i].x)+R121);
+				G211=log((x2-obs_p_[i].x)+R211);
+				G111=log((x1-obs_p_[i].x)+R111);
+
+				temp_obs[j] = G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
+			}
+		}
+
+		for (int n = 0; n < model_num_; n++){
+			obs_p_[i].val += temp_obs[n];
+		}
+	}
+
+	for (int i = 0; i < obs_num_; i++){
+		temp_noise = dist(generator);
+		obs_p_[i].val += temp_noise;
+		obs_p_[i].dev = noise_dev;
+	}
+	return 0;
+}

src/forward_vzz.cpp

@@ -0,0 +1,72 @@
+#include "gm3d.h"
+
+int GM3D::ForwardVzz(char* noise_level){
+	int i,j;
+	double x1,x2,y1,y2,z1,z2;
+	double R222,R122,R212,R112,R221,R121,R211,R111;
+	double G222,G122,G212,G112,G221,G121,G211,G111;
+
+	double noise_mean, noise_dev, temp_noise;
+	if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
+		noise_mean = noise_dev = 0.0;
+	}
+
+	//添加高斯噪声值
+	default_random_engine generator;
+	normal_distribution<double> dist(noise_mean, noise_dev);
+
+	_1dArray temp_obs(model_num_);
+
+	ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
+	for (i = 0; i < obs_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,x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
+		for (j = 0; j < model_num_; j++){
+			if (fabs(forward_model_[j]) > ZERO){
+				x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
+				y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
+				z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
+
+				R222=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R122=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R212=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R112=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R221=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R121=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				G222=log((x2-obs_p_[i].x)+R222);
+				G122=log((x1-obs_p_[i].x)+R122);
+				G212=log((x2-obs_p_[i].x)+R212);
+				G112=log((x1-obs_p_[i].x)+R112);
+				G221=log((x2-obs_p_[i].x)+R221);
+				G121=log((x1-obs_p_[i].x)+R121);
+				G211=log((x2-obs_p_[i].x)+R211);
+				G111=log((x1-obs_p_[i].x)+R111);
+
+				G222=atan((x2-obs_p_[i].x)*(y2-obs_p_[i].y)/(R222*(z2-obs_p_[i].z)));
+				G122=atan((x1-obs_p_[i].x)*(y2-obs_p_[i].y)/(R122*(z2-obs_p_[i].z)));
+				G212=atan((x2-obs_p_[i].x)*(y1-obs_p_[i].y)/(R212*(z2-obs_p_[i].z)));
+				G112=atan((x1-obs_p_[i].x)*(y1-obs_p_[i].y)/(R112*(z2-obs_p_[i].z)));
+				G221=atan((x2-obs_p_[i].x)*(y2-obs_p_[i].y)/(R221*(z1-obs_p_[i].z)));
+				G121=atan((x1-obs_p_[i].x)*(y2-obs_p_[i].y)/(R121*(z1-obs_p_[i].z)));
+				G211=atan((x2-obs_p_[i].x)*(y1-obs_p_[i].y)/(R211*(z1-obs_p_[i].z)));
+				G111=atan((x1-obs_p_[i].x)*(y1-obs_p_[i].y)/(R111*(z1-obs_p_[i].z)));
+
+				temp_obs[j] = -1.0*G0*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
+			}
+		}
+
+		for (int n = 0; n < model_num_; n++){
+			obs_p_[i].val += temp_obs[n];
+		}
+	}
+
+	for (int i = 0; i < obs_num_; i++){
+		temp_noise = dist(generator);
+		obs_p_[i].val += temp_noise;
+		obs_p_[i].dev = noise_dev;
+	}
+	return 0;
+}

src/forward_za.cpp

@@ -0,0 +1,103 @@
+#include "gm3d.h"
+
+int GM3D::ForwardZa(char* noise_level,char* mag_para){
+	int i,j;
+	double I0,A0,I,A;
+	double Alpha,Beta,Gamma;
+	double x1,x2,y1,y2,z1,z2;
+	double w111_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1;
+	double w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2;
+	double R222,R122,R212,R112,R221,R121,R211,R111;
+	double G222,G122,G212,G112,G221,G121,G211,G111;
+
+	//确定噪声水平
+	double noise_mean, noise_dev, temp_noise;
+	if (2 != sscanf(noise_level,"%lf/%lf",&noise_mean,&noise_dev)){
+		noise_mean = noise_dev = 0.0;
+	}
+
+	//确定磁化参数
+	if (4 != sscanf(mag_para,"%lf/%lf/%lf/%lf",&I0,&A0,&I,&A)){
+		I = 90; A = 0;
+	}
+	else{
+		I=I*Pi/180;
+		A=A*Pi/180;
+
+		Alpha=cos(I)*cos(A);
+		Beta=cos(I)*sin(A);
+		Gamma=sin(I);
+	}
+
+	//添加高斯噪声值
+	default_random_engine generator;
+	normal_distribution<double> dist(noise_mean, noise_dev);
+
+	_1dArray temp_obs(model_num_);
+
+	ProgressBar *bar = new ProgressBar(obs_num_,"Forward modeling");
+	for (i = 0; i < obs_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_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1,\
+		w111_2,w112_2,w121_2,w122_2,w211_2,w212_2,w221_2,w222_2,\
+		x1,x2,y1,y2,z1,z2) shared(i) schedule(guided)
+		for (j = 0; j < model_num_; j++){
+			if (fabs(forward_model_[j]) > ZERO){
+				x1 = model_cube_[j].cen.x - 0.5*model_cube_[j].dx; x2 = model_cube_[j].cen.x + 0.5*model_cube_[j].dx;
+				y1 = model_cube_[j].cen.y - 0.5*model_cube_[j].dy; y2 = model_cube_[j].cen.y + 0.5*model_cube_[j].dy;
+				z1 = model_cube_[j].cen.z - 0.5*model_cube_[j].dz; z2 = model_cube_[j].cen.z + 0.5*model_cube_[j].dz;
+
+				R222=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R122=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R212=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R112=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z2-obs_p_[i].z)*(z2-obs_p_[i].z));
+				R221=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R121=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y2-obs_p_[i].y)*(y2-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R211=sqrt((x2-obs_p_[i].x)*(x2-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+				R111=sqrt((x1-obs_p_[i].x)*(x1-obs_p_[i].x)+(y1-obs_p_[i].y)*(y1-obs_p_[i].y)+(z1-obs_p_[i].z)*(z1-obs_p_[i].z));
+
+				w222_1=(x2-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w122_1=(x1-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w212_1=(x2-obs_p_[i].x)*(y1-obs_p_[i].y);
+				w112_1=(x1-obs_p_[i].x)*(y1-obs_p_[i].y);
+				w221_1=(x2-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w121_1=(x1-obs_p_[i].x)*(y2-obs_p_[i].y);
+				w211_1=(x2-obs_p_[i].x)*(y1-obs_p_[i].y);
+				w111_1=(x1-obs_p_[i].x)*(y1-obs_p_[i].y);
+
+				w222_2=R222*(z2-obs_p_[i].z);
+				w122_2=R122*(z2-obs_p_[i].z);
+				w212_2=R212*(z2-obs_p_[i].z);
+				w112_2=R112*(z2-obs_p_[i].z);
+				w221_2=R221*(z1-obs_p_[i].z);
+				w121_2=R121*(z1-obs_p_[i].z);
+				w211_2=R211*(z1-obs_p_[i].z);
+				w111_2=R111*(z1-obs_p_[i].z);
+				
+				G222=-Gamma*atan2(w222_1,w222_2)+Alpha*log(R222+y2-obs_p_[i].y)+Beta*log(R222+x2-obs_p_[i].x);
+				G122=-Gamma*atan2(w122_1,w122_2)+Alpha*log(R122+y2-obs_p_[i].y)+Beta*log(R122+x1-obs_p_[i].x);
+				G212=-Gamma*atan2(w212_1,w212_2)+Alpha*log(R212+y1-obs_p_[i].y)+Beta*log(R212+x2-obs_p_[i].x);
+				G112=-Gamma*atan2(w112_1,w112_2)+Alpha*log(R112+y1-obs_p_[i].y)+Beta*log(R112+x1-obs_p_[i].x);
+				G221=-Gamma*atan2(w221_1,w221_2)+Alpha*log(R221+y2-obs_p_[i].y)+Beta*log(R221+x2-obs_p_[i].x);
+				G121=-Gamma*atan2(w121_1,w121_2)+Alpha*log(R121+y2-obs_p_[i].y)+Beta*log(R121+x1-obs_p_[i].x);
+				G211=-Gamma*atan2(w211_1,w211_2)+Alpha*log(R211+y1-obs_p_[i].y)+Beta*log(R211+x2-obs_p_[i].x);
+				G111=-Gamma*atan2(w111_1,w111_2)+Alpha*log(R111+y1-obs_p_[i].y)+Beta*log(R111+x1-obs_p_[i].x);
+
+				temp_obs[j] = T0/(4*Pi)*(G222-G122-G212+G112-G221+G121+G211-G111)*forward_model_[j];
+			}
+		}
+
+		for (int n = 0; n < model_num_; n++){
+			obs_p_[i].val += temp_obs[n];
+		}
+	}
+
+	for (int i = 0; i < obs_num_; i++){
+		temp_noise = dist(generator);
+		obs_p_[i].val += temp_noise;
+		obs_p_[i].dev = noise_dev;
+	}
+	return 0;
+}

src/gm3d.h

@@ -0,0 +1,58 @@
+#ifndef _GM3D_H
+#define _GM3D_H
+#include "head_func.h"
+#include "progress_bar.h"
+
+class GM3D{
+public:
+	GM3D(){}
+	~GM3D(){}
+	int BuildRegularGrid(char*); //初始化反演模型空间
+	int AddModels(char*); //读取模型块体参数文件
+	int AddRegularBlock(modelist); //添加普通模型块体
+	int AddTiltedBlock(modelist); //添加倾斜模型块体
+	int AddSphereBlock(modelist); //添加球体椭球体块体
+	int AddInterfaceBlock(modelist); //添加密度界面
+	//模型操作
+	int ReadModel(char*,char*);
+	//输出模型
+	int RegisteredOuput(bool); //注册输出的块体模型
+	int OutMshFile(char*,string); //输出模型文件
+	int OutNeighborFile(char*,char*); //输出模型块体或顶点的相邻关系 暂缓
+	//观测数据
+	int InitObs(char*);
+	int OutObs(char*);
+	//正演函数
+	int ForwardVz(char*);
+	int ForwardVzx(char*);
+	int ForwardVzy(char*);
+	int ForwardVzz(char*);
+	int ForwardDeltaT(char*,char*);
+	int ForwardDeltaTx(char*,char*);
+	int ForwardDeltaTy(char*,char*);
+	int ForwardDeltaTz(char*,char*);
+	int ForwardHax(char*,char*);
+	int ForwardHay(char*,char*);
+	int ForwardZa(char*,char*);
+private:
+	int obs_num_, model_num_, vert_num_;
+	//正演数组
+	obspointArray obs_p_;
+	_2dArray input_models_;
+	_1sArray input_model_names_;
+	_1dArray forward_model_;
+	//模型数据
+	cubeArray model_cube_;
+	cpointArray model_vert_;
+	_1dArray model_block_val_;
+	modelistArray model_list_;
+	_1iArray out_ele_ids_;
+	_1iArray out_ele_data_ids_;
+	_1iArray out_vert_ids_;
+	_int2intMap vert_out_map_;
+	_int2intMap ele_data_out_map_;
+
+	_2iArray model_vert_neighbor_; //暂缓
+	_2iArray model_cube_neighbor_; //暂缓
+};
+#endif

src/head_func.cpp

@@ -0,0 +1,95 @@
+#include "head_func.h"
+/*************************数据结构函数********************************/
+//cpoint减法
+cpoint operator -(cpoint a, cpoint b){
+	cpoint m;
+	m.x=a.x-b.x;
+	m.y=a.y-b.y;
+	m.z=a.z-b.z;
+	return m;
+}
+//cpoint模长
+double modCpoint(cpoint v){
+	return sqrt(v.x*v.x+v.y*v.y+v.z*v.z);
+}
+/*************************全局函数********************************/
+//正负分离的atan函数 正数返回atan 负数返回atan+pi
+double arctg(double v){
+	double ang;
+	if(v>=0) ang=atan(v);
+	else if(v<0) ang=atan(v)+Pi;
+	return ang;
+}
+//将string转换为stringstream
+stringstream str2ss(string s){
+	stringstream sstr;
+	sstr.str(""); sstr.clear(); sstr.str(s);
+	return sstr;
+}
+//返回一个cpoint的位置id字符串
+string cpoint_id(cpoint c){
+	string vert_id, mid_id;
+	stringstream sstemp;
+	sstemp.str(""); sstemp.clear(); sstemp<<setprecision(PRECISION)<<c.x;
+	sstemp>>vert_id;
+	sstemp.str(""); sstemp.clear(); sstemp<<setprecision(PRECISION)<<c.y;
+	sstemp>>mid_id;
+	vert_id = vert_id + " " + mid_id;
+	sstemp.str(""); sstemp.clear(); sstemp<<setprecision(PRECISION)<<c.z;
+	sstemp>>mid_id;
+	vert_id = vert_id + " " + mid_id;
+	return vert_id;
+}
+//测试打开输入文件 如果成功则返回0并输出信息 否则返回1
+int open_infile(ifstream &infile,char* filename){
+	infile.open(filename);
+	if (!infile){
+		cerr << BOLDRED << "error ==> " << RESET << "file not found: " << filename << endl;
+		return -1;
+	}
+	return 0;
+}
+//测试打开输出文件 如果成功则返回0并输出信息 否则返回1
+int open_outfile(ofstream &outfile,char* filename){
+	outfile.open(filename);
+	if (!outfile){
+		cerr << BOLDRED << "error ==> " << RESET << "fail to create the file: " << filename << endl;
+		return -1;
+	}
+	return 0;
+}
+
+//规则网络插值 长方形内数据插值 距离平方反比
+/*长方体示意图*/
+// y
+// |
+// |
+// 3------------2
+// |            |
+// |            |
+// |            |
+// 0------------1--->x
+// 左下角坐标x0 y0
+// 块体尺寸dx dy
+// 插值点坐标x y
+// 四个角点值
+double grid_interpolate(double x0,double y0,double dx,double dy,double x,double y,
+					double d0,double d1,double d2,double d3)
+{
+	double res = 0;
+	double total_dist = 0;
+	double dist[4] = {0,0,0,0};
+	double val[4];
+	val[0] = d0; val[1] = d1; val[2] = d2; val[3] = d3;
+	dist[0] = 1.0/(1e-30+(x-x0)*(x-x0)+(y-y0)*(y-y0));
+	dist[1] = 1.0/(1e-30+(x-dx-x0)*(x-dx-x0)+(y-y0)*(y-y0));
+	dist[2] = 1.0/(1e-30+(x-dx-x0)*(x-dx-x0)+(y-dy-y0)*(y-dy-y0));
+	dist[3] = 1.0/(1e-30+(x-x0)*(x-x0)+(y-dy-y0)*(y-dy-y0));
+	for (int i = 0; i < 4; i++){
+		total_dist += dist[i];
+	}
+	for (int i = 0; i < 4; i++){
+		res += val[i]*dist[i]/total_dist;
+	}
+	return res;
+}

src/head_func.h

@@ -0,0 +1,93 @@
+#ifndef _HEAD_FUNC_H
+#define _HEAD_FUNC_H
+#include "iostream"
+#include "fstream"
+#include "sstream"
+#include "string.h"
+#include "cmath"
+#include "iomanip"
+#include "stdio.h"
+#include "stdlib.h"
+#include "unistd.h"
+#include "vector"
+#include "map"
+#include "algorithm"
+#include "ctime"
+#include "omp.h"
+#include "random"
+
+using namespace std;
+
+//数学常量
+#define BDL_MAX 1e+30 ///< 定义变量的最大值
+#define BDL_MIN -1e+30 ///< 定义变量的最小值
+#define PRECISION 16 ///< 定义小数点后需要使用的位数
+#define ZERO 1e-16 ///< 定义零值
+
+//物理常量
+#define Pi (4.0*atan(1.0)) ///< 圆周率
+#define G0 6.67408e-3 ///< 万有引力常数。注意这里的量级本来应该是e-11，考虑到单位转换，取维度单位为m，密度单位为g/cm^3，乘以G0则重力单位即为mGal
+#define T0 5.0e+4 ///< 地磁场平均强度
+//宏函数
+#define MAX(a,b) (a>b?a:b) ///< 返回a与b的最大值
+#define MIN(a,b) (a<b?a:b) ///< 返回a与b的最小值
+#define SetToBox(a,b,in) (MAX(a,MIN(b,in))) ///< 返回a与b之间的一个值，若in在a与b之间则直接返回，否则返回较近的边界值
+//终端显示控制符
+#define BOLDRED "\033[1m\033[31m" ///< 设置后续字符字体为红色加粗
+#define BOLDGREEN "\033[1m\033[32m" ///< 设置后续字符字体为绿色加粗
+#define BOLDYELLOW "\033[1m\033[33m" ///< 设置后续字符字体为黄色加粗
+#define BOLDBLUE "\033[1m\033[34m" ///< 设置后续字符字体为蓝色加粗
+#define UNDERLINE "\033[1m\033[4m" ///< 设置后续字符为添加下划线
+#define RESET "\033[0m" ///< 重置字符设置
+#define MOVEUP(x) printf("\033[%dA", (x)) ///< 将光标向上挪x行
+#define MOVEDOWN(x) printf("\033[%dB", (x)) ///< 将光标向下娜x行
+#define MOVELEFT(x) printf("\033[%dD", (x)) ///< 将光标向左娜x字符
+#define MOVERIGHT(x) printf("\033[%dC", (x)) ///< 将光标向右娜x字符
+#define MOVETO(y,x) printf("\033[%d;%dH", (y), (x)) ///< 将光标向右娜动y字符,向上挪动x字符
+#define CLEARLINE "\033[K" ///< 清除本行
+#define CLEARALL "\033[2J" ///< 清除终端满屏
+//数据结构
+typedef vector<int> _1iArray; ///< 整形一维向量
+typedef vector<double> _1dArray; ///< 双精度浮点一维向量
+typedef vector<string> _1sArray; ///< 字符串一维向量
+typedef vector<vector<int> > _2iArray; ///< 整形浮点二维向量
+typedef vector<vector<double> > _2dArray; ///< 双精度浮点二维向量
+typedef map<int,int> _int2intMap; ///< 整型到整形的映射
+//模型块体参数
+struct modelist{
+	char mod_type[1024];
+	char val_type[1024];
+	char mod_para[1024];
+	double mod_value;
+};
+typedef vector<modelist> modelistArray;
+//直角坐标系点
+struct cpoint{
+	int id = -1;
+	double x = BDL_MAX; double y = BDL_MAX; double z = BDL_MAX;
+};
+typedef vector<cpoint> cpointArray;
+typedef map<string,cpoint> _str2pointMap;
+//观测点
+struct obspoint : public cpoint{
+	double val = BDL_MAX; double dev = BDL_MAX;
+};
+typedef vector<obspoint> obspointArray;
+//块体
+struct cube{
+	cpoint cen;
+	int ids[8] = {-1,-1,-1,-1,-1,-1,-1,-1};
+	double dx = BDL_MAX; double dy = BDL_MAX; double dz = BDL_MAX;
+};
+typedef vector<cube> cubeArray;
+/*************************数据结构函数********************************/
+cpoint operator -(cpoint,cpoint); //矢量减法
+double modCpoint(cpoint); //矢量模
+/*************************全局函数********************************/
+double arctg(double); //正负分离的atan函数 正数返回atan 负数返回atan+pi
+stringstream str2ss(string); //将string转换为stringstream
+string cpoint_id(cpoint); //返回一个cpoint的位置id
+int open_infile(ifstream&,char*); //测试打开输入文件 如果成功则返回0并输出信息 否则返回1
+int open_outfile(ofstream&,char*); //测试打开输出文件 如果成功则返回0并输出信息 否则返回1
+double grid_interpolate(double,double,double,double,double,double,double,double,double,double); //规则网络插值
+#endif

src/init_obs.cpp

@@ -0,0 +1,61 @@
+#include "gm3d.h"
+
+int GM3D::InitObs(char* obs_para){
+	obspoint temp_obs;
+	string temp_str;
+	stringstream temp_ss;
+	double x,y;
+	double xmin,xmax,ymin,ymax;
+	double xs,xe,ys,ye,eleva,dx,dy;
+
+	//按格式解析参数 初始化观测位置 用于正演计算
+	if (7 == sscanf(obs_para,"%lf/%lf/%lf/%lf/%lf/%lf/%lf",&xs,&dx,&xe,&ys,&dy,&ye,&eleva)){
+		xmin = MIN(xs,xe); xmax = MAX(xs,xe);
+		ymin = MIN(ys,ye); ymax = MAX(ys,ye);
+
+		x = xs;
+		while(x >= xmin && x <= xmax){
+			y = ys;
+			while(y >= ymin && y <= ymax){
+				temp_obs.id = obs_p_.size();
+				temp_obs.x = x; temp_obs.y = y; temp_obs.z = -1.0*eleva;
+				temp_obs.val = temp_obs.dev = 0.0;
+				obs_p_.push_back(temp_obs);
+				y += dy;
+			}
+			x += dx;
+		}
+	}
+	//解析失败 按文件读入 用于反演使用或者正演计算
+	else{
+		ifstream infile;
+		if (open_infile(infile,obs_para)) return -1;
+
+		while(getline(infile,temp_str)){
+			if (*(temp_str.begin()) == '#') continue;
+			else{
+				//按每行3个数据解析 初始化为用于正演的观测点
+				if (3 == sscanf(temp_str.c_str(),"%lf %lf %lf",&temp_obs.y,&temp_obs.x,&temp_obs.z)){
+					temp_obs.z *= -1.0;
+					temp_obs.id = obs_p_.size();
+					temp_obs.val = temp_obs.dev = 0.0;
+					obs_p_.push_back(temp_obs);
+				}
+				else{
+					cerr << BOLDYELLOW << "ignored ==> " << RESET << "wrong input: " << temp_str << endl;
+					continue;
+				}
+			}
+		}
+		infile.close();
+	}
+
+	if (obs_p_.empty()){
+		cerr << BOLDRED << "error ==> " << RESET << "fail to initial observations with the parameter: " << obs_para << endl;
+		return -1;
+	}
+	else{
+		obs_num_ = obs_p_.size();
+	}
+	return 0;
+}

src/main.cpp

@@ -0,0 +1,196 @@
+#include "gm3d.h"
+#include "disp_help.h"
+
+void display_help_info(char* program_name){
+	string ex_name = program_name;
+	string ex_usage = program_name;
+	ex_usage += " {-o<output-msh-file> -m<model-list-file>  \
+	[-d<xs>/<dx>/<xe>/<ys>/<dy>/<ye>/<zs>/<dz>/<ze>] [-r]} \
+	| {-i<input-msh-file> -p<observation-file>|<xs>/<dx>/<xe>/<ys>/<dy>/<ye>/<elevation> \
+	-f<output-data-file> -tVz|Vzx|Vzy|Vzz|DT|DTx|DTy|DTz|Hax|Hay|Za \
+	[-v<I0>/<D0>/<I>/<D>] [-n<noise-mean>/<noise-dev>]} -e<element-data-name> [-h]";
+
+	DispHelp dh;
+	dh.changeLayerOut(0,10);
+	dh.addHeadInfo(ex_name,"1.0","3D forward modeling of gravity and magnetic data under the Cartesian coordinates.","Dr. Yi Zhang (zhangyiss@icloud.com). School of Earth Sciences, Zhejiang University");
+	dh.addUsage(ex_usage);
+	dh.addOption("Filename of the input Gmsh(.msh) model file for forward calculation.","-i");
+	dh.addOption("Filename of the output Gmsh(.msh) model file built with given parameters.","-o");
+	dh.addOption("Filename of the output observation file of gravity or magnetic data.","-f");
+	dh.addOption("3D dimensions of the model space. the suffix 's' means the starting coordinate and 'e' represents the ending coordinate in axial directions. \
+	 'dx', 'dy' and 'dz' are step lengths. The default value is 10/20/990/10/20/990/10/20/490. \
+	 The axial orientation adopted by the program is a right-hand Cartesian system with the z-axis point vertical downward.","-d");
+	dh.addOption("Model file that contains different types of model parameter. See instructions for formats of different model types. Note that any line starts with '#' will be skipped. \
+	Four different types of models are supported by the program. Keyword parameters for different model types are shown as:","-m");
+	dh.addOptionSec("1. regular_block add|replace|erase <physical_value> <xmin>/<xmax>/<ymin>/<ymax>/<zmin>/<zmax>");
+	dh.addOptionSec("2. tilted_block add|replace|erase <physical_value> <xmin_z>/<xmax_z>/<ymin_z>/<ymax_z>/<zmin>/<xmin_Z>/<xmax_Z>/<ymin_Z>/<ymax_Z>/<zmax>");
+	dh.addOptionSec("3. sphere add|replace|erase <physical_value> <xcen>/<ycen>/<zcen>/<xradius>/<yradius>/<zradius>");
+	dh.addOptionSec("4. interface add|replace|erase/top|bot <physical_value> <filename>");
+	dh.addOption("Element data name of the input/output Gmsh(.msh) file. Note that the name must be around by \"\".","-e");
+	dh.addOption("Observation locations. You can either initialize the observation points from parameters or a file. \
+	 Each line of the file contain coordinates y(easting), x(northing) and z(elevation) of an observation point.","-p");
+	dh.addOption("Forward component Vz, Vzx, Vzy or Vzz for gravitational data and \
+		DeltaT, DeltaTx, DeltaTy, DeltaTz, Hax, Hay and Za for magnetic data.","-t");
+	dh.addOption("Inclination and declination of the geomagnetic field and magnetization.","-v");
+	dh.addOption("Add noise to the forward calculated data","-n");
+	dh.addOption("Remove model elements with no data in the output Gmsh(.msh) file.","-r");
+	dh.addOption("Display help information.","-h");
+	dh.show();
+	return;
+}
+
+int main(int argc, char* argv[]){
+	GM3D gm3d_instance;
+	char dimension[1024] = "10/20/990/10/20/990/10/20/490";
+	char modelname[1024] = "NULL";
+	char in_mshname[1024] = "NULL";
+	char out_mshname[1024] = "Untitled.msh";
+	char elename[1024] = "Untitled";
+	char obsername[1024] = "NULL";
+	char out_obserfile[1024] = "NULL";
+	char noise_para[1024] = "NULL";
+	char forward_type[1024] = "Vz";
+	char mag_field_para[1024] = "90/0/90/0";
+	bool build_model = true;
+	bool remove_null = false;
+
+	opterr = 0; //内置参数 若不为0则会在发生遭遇错误时输出一条信息到屏幕
+
+	if (argc == 1){
+		display_help_info(argv[0]);
+		return 0;
+	}
+
+	int curr;
+	/*循环拾取参数 最后一个参数为-1 需要变量的参数后跟一个冒号 可有可无参数跟两个冒号*/
+	while((curr = getopt(argc,argv,"hrd:m:i:o:e:p:f:n:t:v:")) != -1){
+		/*匹配命令*/
+		switch (curr){
+			case 'h': //显示帮助信息
+				display_help_info(argv[0]);
+				return 0;
+			case 'r':
+				remove_null = true;
+				break;
+			case 'd':
+				if (1!=sscanf(optarg,"%s",dimension)){
+					cout << "error ==> wrong format of " << optarg << endl;
+				}
+				break;
+			case 'm':
+				if (1!=sscanf(optarg,"%s",modelname)){
+					cout << "error ==> wrong format of " << optarg << endl;
+				}
+				break;
+			case 'i':
+				if (1!=sscanf(optarg,"%s",in_mshname)){
+					cout << "error ==> wrong format of " << optarg << endl;
+				}
+				else build_model = false;
+				break;
+			case 'o':
+				if (1!=sscanf(optarg,"%s",out_mshname)){
+					cout << "error ==> wrong format of " << optarg << endl;
+				}
+				break;
+			case 'f':
+				if (1!=sscanf(optarg,"%s",out_obserfile)){
+					cout << "error ==> wrong format of " << optarg << endl;
+				}
+				break;
+			case 'n':
+				if (1!=sscanf(optarg,"%s",noise_para)){
+					cout << "error ==> wrong format of " << optarg << endl;
+				}
+				break;
+			case 'e':
+				if (1!=sscanf(optarg,"%s",elename)){
+					cout << "error ==> wrong format of " << optarg << endl;
+				}
+				break;
+			case 'v':
+				if (1!=sscanf(optarg,"%s",mag_field_para)){
+					cout << "error ==> wrong format of " << optarg << endl;
+				}
+				break;
+			case 'p':
+				if (1!=sscanf(optarg,"%s",obsername)){
+					cout << "error ==> wrong format of " << optarg << endl;
+				}
+				break;
+			case 't':
+				if (1!=sscanf(optarg,"%s",forward_type)){
+					cout << "error ==> wrong format of " << optarg << endl;
+				}
+				break;
+			case '?': //处理未定义或错误参数
+				if (optopt == 'd' || optopt == 'm' || optopt == 'o' 
+					|| optopt == 'e' || optopt == 'i' || optopt == 'p' 
+					|| optopt == 'n' || optopt == 't' || optopt == 'v'){
+					fprintf (stderr, "Option -%c requires an argument.\n", optopt);
+					return -1;
+				}
+				else if (isprint(optopt)){
+					fprintf (stderr, "Unknown option `-%c'.\n", optopt);
+					return -1;
+				}
+				else{
+					fprintf (stderr,"Unknown option character `\\x%x'.\n",optopt);
+					return -1;
+				}
+				break;
+			default:
+				abort();
+		}
+	}
+
+	if (build_model){
+		if (gm3d_instance.BuildRegularGrid(dimension)) return 0;
+		if (gm3d_instance.AddModels(modelname)) return 0;
+		gm3d_instance.RegisteredOuput(remove_null);
+		if (gm3d_instance.OutMshFile(out_mshname,elename)) return 0;
+	}
+	else{
+		if (gm3d_instance.ReadModel(in_mshname,elename)) return 0;
+		if (gm3d_instance.InitObs(obsername)) return 0;
+		if (!strcmp(forward_type,"Vz")){
+			gm3d_instance.ForwardVz(noise_para);
+		}
+		else if (!strcmp(forward_type,"Vzx")){
+			gm3d_instance.ForwardVzx(noise_para);
+		}
+		else if (!strcmp(forward_type,"Vzy")){
+			gm3d_instance.ForwardVzy(noise_para);
+		}
+		else if (!strcmp(forward_type,"Vzz")){
+			gm3d_instance.ForwardVzz(noise_para);
+		}
+		else if (!strcmp(forward_type,"DeltaT")){
+			gm3d_instance.ForwardDeltaT(noise_para,mag_field_para);
+		}
+		else if (!strcmp(forward_type,"DeltaTx")){
+			gm3d_instance.ForwardDeltaTx(noise_para,mag_field_para);
+		}
+		else if (!strcmp(forward_type,"DeltaTy")){
+			gm3d_instance.ForwardDeltaTy(noise_para,mag_field_para);
+		}
+		else if (!strcmp(forward_type,"DeltaTz")){
+			gm3d_instance.ForwardDeltaTz(noise_para,mag_field_para);
+		}
+		else if (!strcmp(forward_type,"Hax")){
+			gm3d_instance.ForwardHax(noise_para,mag_field_para);
+		}
+		else if (!strcmp(forward_type,"Hay")){
+			gm3d_instance.ForwardHay(noise_para,mag_field_para);
+		}
+		else if (!strcmp(forward_type,"Za")){
+			gm3d_instance.ForwardZa(noise_para,mag_field_para);
+		}
+		else{
+			cerr << BOLDYELLOW << "error ==> " << RESET << "wrong forward component type: " << forward_type << endl;
+			return 0;
+		}
+		if (gm3d_instance.OutObs(out_obserfile)) return 0;
+	}
+	return 0;
+}

src/out_msh_file.cpp

@@ -0,0 +1,37 @@
+#include "gm3d.h"
+
+int GM3D::OutMshFile(char* filename,string data_name){
+	if (!strcmp(filename,"NULL")) return 0;
+
+	ofstream outfile;
+	if (open_outfile(outfile,filename)) return -1;
+
+	//好啦 我们这里输出的模型类型应该是块体
+	outfile<<"$MeshFormat"<<endl<<"2.2 0 8"<<endl<<"$EndMeshFormat"<<endl<<"$Nodes"<<endl<< out_vert_ids_.size() <<endl;
+	//输出顶点 最后一个不输出
+	for (int i = 0; i < out_vert_ids_.size(); i++){
+		outfile << i << " " 
+		<< setprecision(16) << model_vert_[out_vert_ids_[i]].x << " " 
+		<< model_vert_[out_vert_ids_[i]].y << " " 
+		<< model_vert_[out_vert_ids_[i]].z << endl;
+	}
+	outfile<<"$EndNodes"<<endl;
+	outfile<<"$Elements"<<endl<< out_ele_ids_.size() <<endl;
+	for (int i = 0; i < out_ele_ids_.size(); i++){
+		outfile << i <<" 5 1 0";
+		for (int j = 0; j < 8; j++) outfile << " " << vert_out_map_[model_cube_[out_ele_ids_[i]].ids[j]];
+		outfile << endl;
+	}
+	outfile << "$EndElements"<< endl;
+
+	if (out_ele_data_ids_.size() > 0){
+		outfile<<"$ElementData"<<endl;
+		outfile<<1<<endl<<"\""<<data_name<<"\""<<endl<<1<<endl<<0.0<<endl<<3<<endl<<0<<endl<<1<<endl<< out_ele_data_ids_.size() <<endl;
+		for (int i = 0; i < out_ele_data_ids_.size(); i++){
+			outfile << ele_data_out_map_[i] << " " << setprecision(16) << model_block_val_[out_ele_data_ids_[i]] << endl;
+		}
+		outfile<<"$EndElementData"<< endl;
+	}
+	outfile.close();
+	return 0;
+}

src/out_neighbor_file.cpp

@@ -0,0 +1,52 @@
+#include "gm3d.h"
+
+int GM3D::OutNeighborFile(char* v_name,char* b_name){
+	/*
+	if (strcmp(v_name,"NULL") || strcmp(b_name,"NULL")){
+		//整理块体间的相邻关系 先初始化顶点相邻数组
+		model_vert_neighbor_.resize(vert_num_);
+		for (int i = 0; i < vert_num_; i++)
+			model_vert_neighbor_[i].resize(8,-1);
+
+		//遍历所有块体整理顶点相邻关系
+		for (int i = 0; i < model_num_; i++){
+			for (int j = 0; j < 8; j++){
+				model_vert_neighbor_[model_cube_[i].ids[j]][j] = model_cube_[i].cen.id;
+			}
+		}
+
+		//遍历所有顶点相邻关系 所有共点的块体都被认定为相邻块体
+		model_cube_neighbor_.resize(model_num_);
+
+		//循环顶点相邻列表 都不为-1则相互添加
+		for (int i = 0; i < vert_num_; i++){
+			for (int n = 0; n < 8; n++){
+				if (model_vert_neighbor_[i][n] != -1){
+					for (int k = 1; k < 8; k++){
+						if (model_vert_neighbor_[i][(n+k)%8] != -1){
+							model_cube_neighbor_[model_vert_neighbor_[i][n]].push_back(model_vert_neighbor_[i][(n+k)%8]);
+							model_cube_neighbor_[model_vert_neighbor_[i][(n+k)%8]].push_back(model_vert_neighbor_[i][n]);
+						}
+					}
+				}
+			}
+		}
+
+		vector<int>::iterator pos; //整型向量的迭代器
+		for (int i = 0; i < model_num_; i++){
+			sort(model_cube_neighbor_[i].begin(),model_cube_neighbor_[i].end()); //对顶点序列由小到大排序
+			pos = unique(model_cube_neighbor_[i].begin(),model_cube_neighbor_[i].end()); //获取重复序列开始的位置
+			model_cube_neighbor_[i].erase(pos,model_cube_neighbor_[i].end()); //删除重复点
+		}
+
+		//清理数组
+		for (int i = 0; i < vert_num_; i++){
+			model_vert_neighbor_[i].clear();
+			vector <int>().swap(model_vert_neighbor_[i]);
+		}
+		model_vert_neighbor_.clear();
+		vector < vector<int> >().swap(model_vert_neighbor_);
+	}
+	*/
+	return 0;
+}

src/out_obs.cpp

@@ -0,0 +1,13 @@
+#include "gm3d.h"
+
+int GM3D::OutObs(char* filename){
+	ofstream outfile;
+	if (open_outfile(outfile,filename)) return -1;
+	outfile << "# y(m)-easting x(m)-northing ele(m) obs-val(mGal|Eo) obs-dev(mGal|Eo)" << endl;
+	for (int i = 0; i < obs_num_; i++){
+		outfile << obs_p_[i].y << " " << obs_p_[i].x << " " << -1.0*obs_p_[i].z << " "
+		 << setprecision(16) << obs_p_[i].val << " "  << obs_p_[i].dev << endl;
+	}
+	outfile.close();
+	return 0;
+}

src/progress_bar.cpp

@@ -0,0 +1,114 @@
+//#ifdef _WINDOWS
+//#include <windows.h>
+//#else
+//#include <sys/ioctl.h>
+//#endif
+
+#include "progress_bar.h"
+
+ProgressBar::ProgressBar() {}
+
+ProgressBar::ProgressBar(unsigned long n_, const char* description_, std::ostream& out_){
+    
+    n = n_;
+    frequency_update = n_;
+    description = description_;
+    out = &out_;
+    
+    unit_bar = "\u2588";
+    unit_space = "-";
+    desc_width = std::strlen(description);  // character width of description field
+
+}
+
+void ProgressBar::SetFrequencyUpdate(unsigned long frequency_update_){
+    
+    if(frequency_update_ > n){
+        frequency_update = n;    // prevents crash if freq_updates_ > n_
+    }
+    else{
+        frequency_update = frequency_update_;
+    }
+}
+
+void ProgressBar::SetStyle(const char* unit_bar_, const char* unit_space_){
+    
+    unit_bar = unit_bar_;
+    unit_space = unit_space_;
+}
+
+int ProgressBar::GetConsoleWidth(){
+
+    int width;
+
+    #ifdef _WINDOWS
+        CONSOLE_SCREEN_BUFFER_INFO csbi;
+        GetConsoleScreenBufferInfo(GetStdHandle(STD_OUTPUT_HANDLE), &csbi);
+        width = csbi.srWindow.Right - csbi.srWindow.Left;
+    #else
+        struct winsize win;
+        //注意！当我们使用pipe here-doc等通道获取程序参数时无法正确的获取窗口大小 此时我们将使用预定值
+        if (ioctl(0, TIOCGWINSZ, &win) != -1)
+            width = win.ws_col;
+        else width = 100;
+    #endif
+
+    return width;
+}
+
+int ProgressBar::GetBarLength(){
+
+    // get console width and according adjust the length of the progress bar
+
+    int bar_length = static_cast<int>((GetConsoleWidth() - desc_width - CHARACTER_WIDTH_PERCENTAGE) / 2.);
+
+    return bar_length;
+}
+
+void ProgressBar::ClearBarField(){
+
+    for(int i=0;i<GetConsoleWidth();++i){
+        *out << " ";
+    }
+    *out << "\r" << std::flush;
+}
+
+void ProgressBar::Progressed(unsigned long idx_)
+{
+    try{
+        if(idx_ > n) throw idx_;
+
+        // determines whether to update the progress bar from frequency_update
+        if ((idx_ != n-1) && ((idx_+1) % (n/frequency_update) != 0)) return;
+
+        // calculate the size of the progress bar
+        int bar_size = GetBarLength();
+    
+        // calculate percentage of progress
+        double progress_percent = idx_* TOTAL_PERCENTAGE/(n-1);
+
+        // calculate the percentage value of a unit bar 
+        double percent_per_unit_bar = TOTAL_PERCENTAGE/bar_size;
+
+        // display progress bar
+        *out << " " << description << " |";
+
+        for(int bar_length=0;bar_length<=bar_size-1;++bar_length){
+            if(bar_length*percent_per_unit_bar<progress_percent){
+                *out << unit_bar;
+            }
+            else{
+                *out << unit_space;
+            }
+        }
+
+        if(idx_ == n-1)
+            *out << "|" << std::setw(CHARACTER_WIDTH_PERCENTAGE + 1) << std::setprecision(1) << std::fixed << progress_percent << "%\r" << std::flush << std::endl;
+        else *out << "|" << std::setw(CHARACTER_WIDTH_PERCENTAGE + 1) << std::setprecision(1) << std::fixed << progress_percent << "%\r" << std::flush;
+    }
+    catch(unsigned long e){
+        ClearBarField();
+        std::cerr << "PROGRESS_BAR_EXCEPTION: _idx (" << e << ") went out of bounds, greater than n (" << n << ")." << std::endl << std::flush;
+    }
+
+}

src/progress_bar.h

@@ -0,0 +1,41 @@
+#ifndef _PROGRESS_BAR_
+#define _PROGRESS_BAR_
+#include <sys/ioctl.h>
+#include <iostream>
+#include <iomanip>
+#include <cstring>
+#include <thread>
+#include <chrono>
+
+#define TOTAL_PERCENTAGE 100.0
+#define CHARACTER_WIDTH_PERCENTAGE 4
+
+class ProgressBar
+{
+public: 
+
+    ProgressBar();
+    ProgressBar(unsigned long n_, const char *description_="", std::ostream& out_=std::cerr);
+
+    void SetFrequencyUpdate(unsigned long frequency_update_);
+    void SetStyle(const char* unit_bar_, const char* unit_space_);
+
+    void Progressed(unsigned long idx_);
+
+private:
+	
+    unsigned long n;
+    unsigned int desc_width;
+    unsigned long frequency_update;
+    std::ostream* out;
+		
+    const char *description;
+    const char *unit_bar;
+    const char *unit_space;
+		
+    void ClearBarField();
+    int GetConsoleWidth();
+    int GetBarLength();
+
+};
+#endif

src/read_model.cpp

@@ -0,0 +1,122 @@
+#include "gm3d.h"
+
+int GM3D::ReadModel(char* filename,char* input_forward_model_name){
+	int temp_int,ele_type,attri_num,temp_attri,temp_id;
+	double temp_val;
+	_1dArray temp_model;
+	cpoint temp_vert;
+	cube temp_cu;
+	string temp_str;
+	stringstream temp_ss;
+
+	ifstream mshin;
+	if (open_infile(mshin,filename)) return -1; //检查并打开模型文件
+	while(getline(mshin,temp_str)){
+		//读入模型空间顶点集 msh文件版本为2.2
+		if (temp_str == "$Nodes"){
+			getline(mshin,temp_str);
+			temp_ss = str2ss(temp_str);
+			temp_ss >> vert_num_; //第一个数为顶点的个数
+			model_vert_.resize(vert_num_); //开辟空间
+			for (int i = 0; i < vert_num_; i++){
+				getline(mshin,temp_str);
+				temp_ss = str2ss(temp_str);
+				temp_ss >> temp_vert.id >> temp_vert.x >> temp_vert.y >> temp_vert.z;
+				model_vert_[i] = temp_vert;
+			}
+		}
+		//读入模型空间单元体
+		else if (temp_str == "$Elements"){
+			getline(mshin,temp_str);
+			temp_ss = str2ss(temp_str);
+			temp_ss >> model_num_; //第一个数为总元素的个数 包含了所有类型的元素 比如三角形 四边形 块体等
+			model_cube_.resize(model_num_);
+			for (int i = 0; i < model_num_; i++){
+				getline(mshin,temp_str);
+				temp_ss = str2ss(temp_str);
+				temp_ss >> temp_cu.cen.id >> ele_type;
+				//只读入块体
+				if (ele_type == 5){
+					temp_ss >> attri_num; //跳过模型单元的几何组与物理组等信息 以后可能会有用
+					for (int a = 0; a < attri_num; a++)
+						temp_ss >> temp_attri;
+					for (int a = 0; a < 8; a++)
+						temp_ss >> temp_cu.ids[a];
+					model_cube_[i] = temp_cu;
+				}
+			}
+		}
+		else continue; //不能识别的单元都被忽略了
+	}
+	mshin.close();
+
+	//第二次读入模型文件 初始化模型单元属性
+	if (open_infile(mshin,filename)) return -1; //检查并打开模型文件
+	while(getline(mshin,temp_str)){
+		//读入模型单元属性 注意因为msh文件中$ElementData并未注明所属元素类型
+		//所以可能会将其他元素类型的属性值也读入 但因为其在pyIdMap中并未注册 所以属性值会全为0 在后续使用时我们需要通过名称辨别
+		if (temp_str == "$ElementData"){
+			temp_model.resize(model_num_,0.0); //初始化temp_model 为读入模型单元属性做准备
+			for (int i = 0; i < 2; i++) //先读入元素块的名称 添加到数组
+				getline(mshin,temp_str);
+			input_model_names_.push_back(temp_str);
+			for (int i = 0; i < 6; i++) //跳过元素属性前面的值 最后一次为当前元素块的个数
+				getline(mshin,temp_str);
+			temp_ss = str2ss(temp_str);
+			temp_ss >> temp_int;
+			for (int i = 0; i < temp_int; i++){
+				getline(mshin,temp_str);
+				temp_ss = str2ss(temp_str);
+				temp_ss >> temp_id >> temp_val; //读入单元体索引与属性值
+				temp_model[temp_id] = temp_val;
+			}
+			input_models_.push_back(temp_model);
+			temp_model.clear();
+		}
+		else continue; //不能识别的单元都被忽略了
+	}
+	mshin.close();
+
+	//清理映射
+	temp_model.clear();
+	vector <double>().swap(temp_model);
+
+	//初始化数组
+	forward_model_.resize(model_num_,0.0);
+	//匹配数据名称
+	char forward_model_name[1024] = "\"";
+	strcat(forward_model_name,input_forward_model_name);
+	strcat(forward_model_name,"\"");
+	for (int i = 0; i < input_model_names_.size(); i++){
+		if (!strcmp(input_model_names_[i].c_str(),forward_model_name)){
+			forward_model_ = input_models_[i];
+			clog << "forward model initialized from file." << endl;
+		}
+	}
+
+	//计算块体的中心位置和尺寸
+	cpoint corner[8];
+	for (int i = 0; i < model_num_; i++){
+		for (int j = 0; j < 8; j++){
+			corner[j] = model_vert_[model_cube_[i].ids[j]];
+		}
+
+		model_cube_[i].cen.x = 0.5*(corner[0].x + corner[6].x);
+		model_cube_[i].cen.y = 0.5*(corner[0].y + corner[6].y);
+		model_cube_[i].cen.z = 0.5*(corner[0].z + corner[6].z);
+		model_cube_[i].dx = fabs(corner[6].x - corner[0].x);
+		model_cube_[i].dy = fabs(corner[6].y - corner[0].y);
+		model_cube_[i].dz = fabs(corner[6].z - corner[0].z);
+	}
+
+	for (int i = 0; i < input_models_.size(); i++){
+		input_models_[i].clear();
+		vector <double>().swap(input_models_[i]);
+	}
+	input_models_.clear();
+	vector < vector<double> >().swap(input_models_);
+
+	input_model_names_.clear();
+	vector <string>().swap(input_model_names_);
+	return 0;
+}

src/registered_output.cpp

@@ -0,0 +1,61 @@
+#include "gm3d.h"
+
+int GM3D::RegisteredOuput(bool remove_empty_element){
+	int count;
+	//统计输出模型单元块体和顶点 以及输出的块体数据列表
+	if (remove_empty_element){
+		count = 0;
+		//遍历所有块体数据 注册有值的块体
+		for (int i = 0; i < model_num_; i++){
+			if (model_block_val_[i] != BDL_MAX){
+				out_ele_data_ids_.push_back(i);
+				out_ele_ids_.push_back(i);
+				ele_data_out_map_[count] = count;
+				count++;
+			}
+		}
+
+		//遍历所有注册的块体 添加顶点
+		for (int i = 0; i < out_ele_ids_.size(); i++){
+			for (int j = 0; j < 8; j++){
+				out_vert_ids_.push_back(model_cube_[out_ele_ids_[i]].ids[j]);
+			}
+		}
+
+		//去除输出顶点中的重复部分
+		vector<int>::iterator pos; //整型向量的迭代器
+		sort(out_vert_ids_.begin(),out_vert_ids_.end()); //对顶点序列由小到大排序
+		pos = unique(out_vert_ids_.begin(),out_vert_ids_.end()); //获取重复序列开始的位置
+		out_vert_ids_.erase(pos,out_vert_ids_.end()); //删除重复点
+
+		//将需要输出的模型顶点序号与它们的排序做一个对应 保证在输出文件中顶点索引号始终是从0开始的连续的序列
+		for (int i = 0; i < out_vert_ids_.size(); i++){
+			vert_out_map_[out_vert_ids_[i]] = i;
+		}
+	}
+	else{
+		//输出的模型块体为所有
+		out_ele_ids_.resize(model_num_);
+		for (int i = 0; i < model_num_; i++){
+			out_ele_ids_[i] = i;
+		}
+
+		//输出的模型顶点为所有
+		out_vert_ids_.resize(vert_num_);
+		for (int i = 0; i < vert_num_; i++){
+			vert_out_map_[i] = out_vert_ids_[i] = i;
+		}
+
+		//注册所有有值的块体数据
+		count = 0;
+		for (int i = 0; i < model_num_; i++){
+			if (model_block_val_[i] != BDL_MAX){
+				ele_data_out_map_[count] = i;
+				out_ele_data_ids_.push_back(i);
+				count++;
+			}
+		}
+	}
+	//输出所有模型单元块体和顶点 只统计输出的块体数据列表
+	return 0;
+}