initial update

.gitignore

@@ -30,3 +30,7 @@
 *.exe
 *.out
 *.app
+
+# Mac folder
+.DS_Store
+build_mac/

CMakeLists.txt

@@ -0,0 +1,7 @@
+cmake_minimum_required(VERSION 3.15.2)
+
+set(CMAKE_CXX_COMPILER /usr/local/bin/g++-9)
+
+project(GM3D_GUI VERSION 0.2 LANGUAGES CXX)
+
+add_subdirectory(src)

README.md

@@ -0,0 +1,111 @@
+# gm3d-gui程序说明书
+
+gm3d是一个开源的直角坐标系下的三维重磁场建模与正演的命令行工具。它能够以方块为基本的建模单元构建复杂的三维密度或磁化率模型并正演计算模型的重磁数据及各方向梯度数据等。
+
+## 程序界面说明
+
+### 建模界面
+
+![建模界面](assert/fig-1.png)
+
+![模型添加界面](assert/fig-2.png)
+
+### 正演界面
+
+![正演界面](assert/fig-3.png)
+
+## 参数与文件说明
+
+### Msh本文文件格式
+
+文本格式Gmsh模型文件（.msh）主要格式如下所示。在格式描述中变量名称以尖括号包围，默认值或数据类型以等号给出。格式注释以！！开头（注意实际文件中无此特性，仅为方便说明之用）。
+
+
+```
+！！网格文件类型 文件类型=0表示为ASCII码文本文件 （必需）
+$MeshFormat
+<版本号=2.2> <文件类型=0> <浮点字符长度=8>
+$EndMeshFormat
+！！物理组名称 （非必需）
+$PhysicalNames
+<名称个数=int>
+<物理维度=1|2|3> <物理组标签=int> <物理组名称="string">
+...
+$EndPhysicalNames
+！！点描述（必需）
+$Nodes
+<点个数=int>
+<点编号=int> <x坐标=double> <y坐标=double> <z坐标=double>
+...
+$EndNodes
+！！几何元素描述（必需） 不同几何元素的类型编号请参考Gmsh说明文档 第一个标签为物理组 第二个标签为几何组
+！！元素顶点索引的排序请参考Gmsh说明文档
+$Elements
+<几何元素的个数=int>
+<元素编号=int> <元素类型编码=int> <标签个数=int> <标签=int> <标签=int>... <元素顶点索引=int> <元素顶点索引=int>...
+...
+$EndElements
+！！点数据描述 字符串标签一般即为此数据的名称（非必需）
+$NodeData
+<字符串标签个数=1>
+<字符串标签="string">                           ！！数据名称
+<浮点标签个数=1>
+<浮点标签=double>                              ！！时间戳
+<整形标签个数=3>
+<整形标签=int>                                 ！！时间排序
+<整形标签=1>                                   ！！一维标量
+<整形标签=int>                                 ！！数据个数
+<点索引=int> <数据值=double>
+...
+...
+$EndNodeData
+！！元素数据描述 字符串标签一般即为此数据的名称（非必需）
+$ElementData
+<字符串标签个数=1>
+<字符串标签="string">                           ！！数据名称
+<浮点标签个数=1>
+<浮点标签=double>                              ！！时间戳
+<整形标签个数=3>
+<整形标签=int>                                 ！！时间排序
+<整形标签=1>                                   ！！一维标量
+<整形标签=int>                                 ！！数据个数
+<元素索引=int> <数据值=double>
+...
+...
+$EndElementData
+！！元素顶点数据描述 字符串标签一般即为此数据的名称（非必需）
+$ElementNodeData
+<字符串标签个数=1>
+<字符串标签="string">                           ！！数据名称
+<浮点标签个数=1>
+<浮点标签=double>                              ！！时间戳
+<整形标签个数=3>
+<整形标签=int>                                 ！！时间排序
+<整形标签=1>                                   ！！一维标量
+<整形标签=int>                                 ！！数据个数
+<元素索引=int> <元素的顶点数据个数> <数据值=double> <数据值=double>...
+...
+...
+$EndElementNodeData
+```
+
+### 模型参数文件
+
+模型参数文件中各行均为一个完整的模型体的信息描述，包含以下部分：
+```
+<模型类型> <赋值类型> <物理属性值> <几何属性参数>
+```
+其中：
+
+* <模型类型>包括`regular_block`，`tilted_block`，`sphere`和`interface`四种类型。
+* <赋值类型>包含`add`，`replace`和`erase`三种类型。特别的，对于`interface`类型的模型我们还需要制定操作的对象为界面上`top`或下`bot`，并与前面的赋值类型通过斜杠链接。
+* <物理属性值>为模型体的密度或磁化率值。
+* <几何属性参数>模型体的几何参数依据模型类型各不相同，详见帮助信息。
+
+### 观测点文件
+
+观测点是一个简单的文本文件，文件每行包含观测点的y（东）、x（北）与z（高程）坐标。
+
+### 界面文件
+
+界面文件为一个简单的文本描述的网格文件，文件头需包含"#range=\<xmin\>/\<dx\>/\<xmax\>/\<ymin\>/\<dx\>/\<ymax\>"信息。文件每行包含高程点的y（东）、x（北）与z（高程）坐标。

fluid_project/add_model.cxx

@@ -0,0 +1,94 @@
+// generated by Fast Light User Interface Designer (fluid) version 1.0305
+
+#include "add_model.h"
+
+Fl_Double_Window *add_mod_win=(Fl_Double_Window *)0;
+
+Fl_Group *mod_type_group=(Fl_Group *)0;
+
+Fl_Round_Button *reg_bok_rbtn=(Fl_Round_Button *)0;
+
+Fl_Round_Button *til_bok_rbtn=(Fl_Round_Button *)0;
+
+Fl_Round_Button *sph_rbtn=(Fl_Round_Button *)0;
+
+Fl_Round_Button *int_face_rbtn=(Fl_Round_Button *)0;
+
+Fl_Group *val_type_group=(Fl_Group *)0;
+
+Fl_Round_Button *app_val_rbtn=(Fl_Round_Button *)0;
+
+Fl_Round_Button *rep_val_rbtn=(Fl_Round_Button *)0;
+
+Fl_Round_Button *era_val_rbtn=(Fl_Round_Button *)0;
+
+Fl_Group *agn_part_group=(Fl_Group *)0;
+
+Fl_Round_Button *top_val_btn=(Fl_Round_Button *)0;
+
+Fl_Round_Button *bot_val_btn=(Fl_Round_Button *)0;
+
+Fl_Input *sig_mod_para_input=(Fl_Input *)0;
+
+Fl_Return_Button *can_add_btn=(Fl_Return_Button *)0;
+
+Fl_Button *sig_add_btn=(Fl_Button *)0;
+
+Fl_Input *mod_val_input=(Fl_Input *)0;
+
+void cb_add_mod_btn(Fl_Button*, void*) {
+  { add_mod_win = new Fl_Double_Window(315, 300, "Add model (gm3d)");
+    { mod_type_group = new Fl_Group(20, 25, 250, 53, "Model Type :");
+      mod_type_group->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+      { reg_bok_rbtn = new Fl_Round_Button(20, 25, 110, 28, "Regular Block");
+        reg_bok_rbtn->down_box(FL_ROUND_DOWN_BOX);
+      } // Fl_Round_Button* reg_bok_rbtn
+      { til_bok_rbtn = new Fl_Round_Button(160, 25, 100, 28, "Tilted Block");
+        til_bok_rbtn->down_box(FL_ROUND_DOWN_BOX);
+      } // Fl_Round_Button* til_bok_rbtn
+      { sph_rbtn = new Fl_Round_Button(20, 50, 70, 28, "Sphere");
+        sph_rbtn->down_box(FL_ROUND_DOWN_BOX);
+      } // Fl_Round_Button* sph_rbtn
+      { int_face_rbtn = new Fl_Round_Button(160, 50, 80, 28, "Interface");
+        int_face_rbtn->down_box(FL_ROUND_DOWN_BOX);
+      } // Fl_Round_Button* int_face_rbtn
+      mod_type_group->end();
+    } // Fl_Group* mod_type_group
+    { val_type_group = new Fl_Group(20, 100, 250, 30, "Value Type :");
+      val_type_group->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+      { app_val_rbtn = new Fl_Round_Button(20, 100, 80, 28, "Append");
+        app_val_rbtn->down_box(FL_ROUND_DOWN_BOX);
+      } // Fl_Round_Button* app_val_rbtn
+      { rep_val_rbtn = new Fl_Round_Button(110, 100, 80, 28, "Replace");
+        rep_val_rbtn->down_box(FL_ROUND_DOWN_BOX);
+      } // Fl_Round_Button* rep_val_rbtn
+      { era_val_rbtn = new Fl_Round_Button(200, 100, 60, 28, "Erase");
+        era_val_rbtn->down_box(FL_ROUND_DOWN_BOX);
+      } // Fl_Round_Button* era_val_rbtn
+      val_type_group->end();
+    } // Fl_Group* val_type_group
+    { agn_part_group = new Fl_Group(20, 150, 160, 28, "Assgin Part :");
+      agn_part_group->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+      { top_val_btn = new Fl_Round_Button(20, 150, 50, 28, "Top");
+        top_val_btn->down_box(FL_ROUND_DOWN_BOX);
+        top_val_btn->deactivate();
+      } // Fl_Round_Button* top_val_btn
+      { bot_val_btn = new Fl_Round_Button(110, 150, 70, 28, "Bottom");
+        bot_val_btn->down_box(FL_ROUND_DOWN_BOX);
+        bot_val_btn->deactivate();
+      } // Fl_Round_Button* bot_val_btn
+      agn_part_group->end();
+    } // Fl_Group* agn_part_group
+    { sig_mod_para_input = new Fl_Input(20, 200, 275, 28, "Model Parameter :");
+      sig_mod_para_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+    } // Fl_Input* sig_mod_para_input
+    { can_add_btn = new Fl_Return_Button(210, 255, 85, 28, "Cancel");
+    } // Fl_Return_Button* can_add_btn
+    { sig_add_btn = new Fl_Button(140, 255, 60, 28, "Add");
+    } // Fl_Button* sig_add_btn
+    { mod_val_input = new Fl_Input(20, 255, 110, 28, "Model Value :");
+      mod_val_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+    } // Fl_Input* mod_val_input
+    add_mod_win->end();
+  } // Fl_Double_Window* add_mod_win
+}

fluid_project/add_model.fl

@@ -0,0 +1,79 @@
+# data file for the Fltk User Interface Designer (fluid)
+version 1.0305
+header_name {.h}
+code_name {.cxx}
+Function {cb_add_mod_btn(Fl_Button*, void*)} {open return_type void
+} {
+  Fl_Window add_mod_win {
+    label {Add model (gm3d)} open
+    xywh {735 331 315 300} type Double visible
+  } {
+    Fl_Group mod_type_group {
+      label {Model Type :} open
+      xywh {20 25 250 53} align 5
+    } {
+      Fl_Round_Button reg_bok_rbtn {
+        label {Regular Block}
+        xywh {20 25 110 28} down_box ROUND_DOWN_BOX
+      }
+      Fl_Round_Button til_bok_rbtn {
+        label {Tilted Block}
+        xywh {160 25 100 28} down_box ROUND_DOWN_BOX
+      }
+      Fl_Round_Button sph_rbtn {
+        label Sphere
+        xywh {20 50 70 28} down_box ROUND_DOWN_BOX
+      }
+      Fl_Round_Button int_face_rbtn {
+        label Interface
+        xywh {160 50 80 28} down_box ROUND_DOWN_BOX
+      }
+    }
+    Fl_Group val_type_group {
+      label {Value Type :} open
+      xywh {20 100 250 30} align 5
+    } {
+      Fl_Round_Button app_val_rbtn {
+        label Append
+        xywh {20 100 80 28} down_box ROUND_DOWN_BOX
+      }
+      Fl_Round_Button rep_val_rbtn {
+        label Replace
+        xywh {110 100 80 28} down_box ROUND_DOWN_BOX
+      }
+      Fl_Round_Button era_val_rbtn {
+        label Erase
+        xywh {200 100 60 28} down_box ROUND_DOWN_BOX
+      }
+    }
+    Fl_Group agn_part_group {
+      label {Assgin Part :} open
+      xywh {20 150 160 28} align 5
+    } {
+      Fl_Round_Button top_val_btn {
+        label Top
+        xywh {20 150 50 28} down_box ROUND_DOWN_BOX deactivate
+      }
+      Fl_Round_Button bot_val_btn {
+        label Bottom
+        xywh {110 150 70 28} down_box ROUND_DOWN_BOX deactivate
+      }
+    }
+    Fl_Input sig_mod_para_input {
+      label {Model Parameter :}
+      xywh {20 200 275 28} align 5
+    }
+    Fl_Return_Button can_add_btn {
+      label Cancel
+      xywh {210 255 85 28}
+    }
+    Fl_Button sig_add_btn {
+      label Add
+      xywh {140 255 60 28}
+    }
+    Fl_Input mod_val_input {
+      label {Model Value :}
+      xywh {20 255 110 28} align 5
+    }
+  }
+}

fluid_project/add_model.h

@@ -0,0 +1,30 @@
+// generated by Fast Light User Interface Designer (fluid) version 1.0305
+
+#ifndef add_model_h
+#define add_model_h
+#include <FL/Fl.H>
+#include <FL/Fl_Double_Window.H>
+extern Fl_Double_Window *add_mod_win;
+#include <FL/Fl_Group.H>
+extern Fl_Group *mod_type_group;
+#include <FL/Fl_Round_Button.H>
+extern Fl_Round_Button *reg_bok_rbtn;
+extern Fl_Round_Button *til_bok_rbtn;
+extern Fl_Round_Button *sph_rbtn;
+extern Fl_Round_Button *int_face_rbtn;
+extern Fl_Group *val_type_group;
+extern Fl_Round_Button *app_val_rbtn;
+extern Fl_Round_Button *rep_val_rbtn;
+extern Fl_Round_Button *era_val_rbtn;
+extern Fl_Group *agn_part_group;
+extern Fl_Round_Button *top_val_btn;
+extern Fl_Round_Button *bot_val_btn;
+#include <FL/Fl_Input.H>
+extern Fl_Input *sig_mod_para_input;
+#include <FL/Fl_Return_Button.H>
+extern Fl_Return_Button *can_add_btn;
+#include <FL/Fl_Button.H>
+extern Fl_Button *sig_add_btn;
+extern Fl_Input *mod_val_input;
+void cb_add_mod_btn(Fl_Button*, void*);
+#endif

fluid_project/gm3d_gui.cpp

@@ -0,0 +1,69 @@
+#include "gm3d_gui.h"
+
+Fl_Double_Window *main_window=(Fl_Double_Window *)0;
+Fl_Tabs *main_tabs=(Fl_Tabs *)0;
+Fl_Group *model_tab=(Fl_Group *)0;
+Fl_Input *mesh_para_input=(Fl_Input *)0;
+Fl_Button *mesh_file_btn=(Fl_Button *)0;
+Fl_Button *mod_para_file_btn=(Fl_Button *)0;
+Fl_Input *mod_ele_input_build=(Fl_Input *)0;
+Fl_Button *build_mod_btn=(Fl_Button *)0;
+Fl_Output *mesh_para_output=(Fl_Output *)0;
+Fl_Button *add_mod_btn=(Fl_Button *)0;
+Fl_Button *del_mod_btn=(Fl_Button *)0;
+Fl_Check_Button *rm_emp_bok_check=(Fl_Check_Button *)0;
+Fl_Browser *mod_para_brw=(Fl_Browser *)0;
+Fl_Button *mod_file_out_btn=(Fl_Button *)0;
+Fl_Output *mod_out_file_output=(Fl_Output *)0;
+Fl_Group *forward_tab=(Fl_Group *)0;
+Fl_Group *grav_group=(Fl_Group *)0;
+Fl_Check_Button *Vz_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Vzx_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Vzy_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Vzz_check=(Fl_Check_Button *)0;
+Fl_Button *mod_file_btn=(Fl_Button *)0;
+Fl_Button *obs_file_btn=(Fl_Button *)0;
+Fl_Input *mod_file_input=(Fl_Input *)0;
+Fl_Output *mod_file_output=(Fl_Output *)0;
+Fl_Input *obs_file_input=(Fl_Input *)0;
+Fl_Output *obs_file_output=(Fl_Output *)0;
+Fl_Input *mod_ele_input=(Fl_Input *)0;
+Fl_Group *mag_group=(Fl_Group *)0;
+Fl_Check_Button *DeltaT_check=(Fl_Check_Button *)0;
+Fl_Check_Button *DeltaTx_check=(Fl_Check_Button *)0;
+Fl_Check_Button *DeltaTy_check=(Fl_Check_Button *)0;
+Fl_Check_Button *DeltaTz_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Hax_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Hay_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Za_check=(Fl_Check_Button *)0;
+Fl_Input *mag_para_input=(Fl_Input *)0;
+Fl_Check_Button *mag_data_check=(Fl_Check_Button *)0;
+Fl_Button *cal_btn=(Fl_Button *)0;
+Fl_Input *noise_para_input=(Fl_Input *)0;
+Fl_Check_Button *noise_check=(Fl_Check_Button *)0;
+Fl_Input *res_file_input=(Fl_Input *)0;
+Fl_Button *res_file_btn=(Fl_Button *)0;
+Fl_Output *res_out_file_output=(Fl_Output *)0;
+
+void cb_mesh_para_input(Fl_Input*, void*){}
+void cb_mesh_file_btn(Fl_Widget*, void*){}
+void cb_add_mod_btn(Fl_Widget*, void*){}
+void cb_mod_para_file_btn(Fl_Widget*, void*){}
+void cb_del_mod_btn(Fl_Widget*, void*){}
+void cb_mod_para_brw(Fl_Browser*, void*){}
+void cb_mod_ele_input(Fl_Input*, void*){}
+void cb_rm_emp_bok_check(Fl_Check_Button*, void*){}
+void cb_mod_file_out_btn(Fl_Widget*, void*){}
+void cb_build_mod_btn(Fl_Widget*, void*){}
+void cb_mod_file_input(Fl_Input*, void*){}
+void cb_mod_file_btn(Fl_Widget*, void*){}
+void cb_obs_file_input(Fl_Input*, void*){}
+void cb_obs_file_btn(Fl_Widget*, void*){}
+void cb_res_file_input(Fl_Input*, void*){}
+void cb_res_file_btn(Fl_Widget*, void*){}
+void cb_nosie_check(Fl_Check_Button*, void*){}
+void cb_noise_para_input(Fl_Input*, void*){}
+void cb_mag_data_check(Fl_Check_Button*, void*){}
+void cb_mag_para_input(Fl_Input*, void*){}
+void cb_cal_btn(Fl_Button*, void*){}
+void cb_noise_check(Fl_Check_Button*, void*){}

fluid_project/gm3d_gui.cxx

@@ -0,0 +1,169 @@
+// generated by Fast Light User Interface Designer (fluid) version 1.0305
+
+#include "gm3d_gui.h"
+
+int main(int argc, char **argv) {
+  { main_window = new Fl_Double_Window(500, 600, "gm3d");
+    { main_tabs = new Fl_Tabs(10, 10, 480, 580);
+      { model_tab = new Fl_Group(10, 40, 480, 550, "Build Model");
+        model_tab->hide();
+        { mesh_para_input = new Fl_Input(40, 70, 300, 28, "Input Mesh Parameters :");
+          mesh_para_input->callback((Fl_Callback*)cb_mesh_para_input);
+          mesh_para_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Input* mesh_para_input
+        { mesh_file_btn = new Fl_Button(360, 70, 100, 28, "Mesh File");
+          mesh_file_btn->callback((Fl_Callback*)cb_mesh_file_btn);
+        } // Fl_Button* mesh_file_btn
+        { mod_para_file_btn = new Fl_Button(170, 160, 160, 28, "Add Model From File");
+          mod_para_file_btn->callback((Fl_Callback*)cb_mod_para_file_btn);
+        } // Fl_Button* mod_para_file_btn
+        { mod_ele_input_build = new Fl_Input(40, 445, 215, 28, "Input Model Element Data Name:");
+          mod_ele_input_build->callback((Fl_Callback*)cb_mod_ele_input);
+          mod_ele_input_build->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Input* mod_ele_input_build
+        { build_mod_btn = new Fl_Button(40, 545, 420, 28, "Build Model !");
+          build_mod_btn->callback((Fl_Callback*)cb_build_mod_btn);
+        } // Fl_Button* build_mod_btn
+        { mesh_para_output = new Fl_Output(40, 120, 300, 28, "Mesh Parameters :");
+          mesh_para_output->box(FL_NO_BOX);
+          mesh_para_output->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Output* mesh_para_output
+        { add_mod_btn = new Fl_Button(40, 160, 100, 28, "Add Model");
+          add_mod_btn->callback((Fl_Callback*)cb_add_mod_btn);
+        } // Fl_Button* add_mod_btn
+        { del_mod_btn = new Fl_Button(360, 160, 100, 28, "Delete Model");
+          del_mod_btn->callback((Fl_Callback*)cb_del_mod_btn);
+          del_mod_btn->deactivate();
+        } // Fl_Button* del_mod_btn
+        { rm_emp_bok_check = new Fl_Check_Button(280, 447, 180, 28, "Remove Empty Blocks");
+          rm_emp_bok_check->down_box(FL_DOWN_BOX);
+          rm_emp_bok_check->callback((Fl_Callback*)cb_rm_emp_bok_check);
+        } // Fl_Check_Button* rm_emp_bok_check
+        { mod_para_brw = new Fl_Browser(40, 205, 420, 210);
+          mod_para_brw->callback((Fl_Callback*)cb_mod_para_brw);
+        } // Fl_Browser* mod_para_brw
+        { mod_file_out_btn = new Fl_Button(360, 500, 100, 28, "Model File");
+          mod_file_out_btn->callback((Fl_Callback*)cb_mod_file_out_btn);
+        } // Fl_Button* mod_file_out_btn
+        { mod_out_file_output = new Fl_Output(40, 500, 300, 28, "Output File Name :");
+          mod_out_file_output->box(FL_NO_BOX);
+          mod_out_file_output->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Output* mod_out_file_output
+        model_tab->end();
+      } // Fl_Group* model_tab
+      { forward_tab = new Fl_Group(10, 40, 480, 550, "Forward Modeling");
+        { grav_group = new Fl_Group(40, 340, 190, 58, "Forward Gravitational Data :");
+          grav_group->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+          { Vz_check = new Fl_Check_Button(40, 340, 50, 28, "Vz");
+            Vz_check->down_box(FL_DOWN_BOX);
+          } // Fl_Check_Button* Vz_check
+          { Vzx_check = new Fl_Check_Button(105, 340, 50, 28, "Vzx");
+            Vzx_check->down_box(FL_DOWN_BOX);
+          } // Fl_Check_Button* Vzx_check
+          { Vzy_check = new Fl_Check_Button(170, 340, 50, 28, "Vzy");
+            Vzy_check->down_box(FL_DOWN_BOX);
+          } // Fl_Check_Button* Vzy_check
+          { Vzz_check = new Fl_Check_Button(40, 370, 50, 28, "Vzz");
+            Vzz_check->down_box(FL_DOWN_BOX);
+          } // Fl_Check_Button* Vzz_check
+          grav_group->end();
+        } // Fl_Group* grav_group
+        { mod_file_btn = new Fl_Button(360, 70, 100, 28, "Model File");
+          mod_file_btn->callback((Fl_Callback*)cb_mod_file_btn);
+        } // Fl_Button* mod_file_btn
+        { obs_file_btn = new Fl_Button(360, 130, 100, 28, "Observe File");
+          obs_file_btn->callback((Fl_Callback*)cb_obs_file_btn);
+        } // Fl_Button* obs_file_btn
+        { mod_file_input = new Fl_Input(40, 70, 300, 28, "Input Model FIle :");
+          mod_file_input->callback((Fl_Callback*)cb_mod_file_input);
+          mod_file_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Input* mod_file_input
+        { mod_file_output = new Fl_Output(40, 240, 200, 28, "Chosen Model File :");
+          mod_file_output->box(FL_NO_BOX);
+          mod_file_output->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Output* mod_file_output
+        { obs_file_input = new Fl_Input(40, 130, 300, 28, "Input Observe  FIle :");
+          obs_file_input->callback((Fl_Callback*)cb_obs_file_input);
+          obs_file_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Input* obs_file_input
+        { obs_file_output = new Fl_Output(260, 240, 200, 28, "Chosen Observe File :");
+          obs_file_output->box(FL_NO_BOX);
+          obs_file_output->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Output* obs_file_output
+        { mod_ele_input = new Fl_Input(260, 315, 200, 28, "Model Element Data Name :");
+          mod_ele_input->callback((Fl_Callback*)cb_mod_ele_input);
+          mod_ele_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Input* mod_ele_input
+        { mag_group = new Fl_Group(40, 440, 420, 90);
+          { DeltaT_check = new Fl_Check_Button(40, 475, 65, 25, "DeltaT");
+            DeltaT_check->down_box(FL_DOWN_BOX);
+            DeltaT_check->deactivate();
+          } // Fl_Check_Button* DeltaT_check
+          { DeltaTx_check = new Fl_Check_Button(120, 475, 70, 25, "DeltaTx");
+            DeltaTx_check->down_box(FL_DOWN_BOX);
+            DeltaTx_check->deactivate();
+          } // Fl_Check_Button* DeltaTx_check
+          { DeltaTy_check = new Fl_Check_Button(200, 475, 70, 25, "DeltaTy");
+            DeltaTy_check->down_box(FL_DOWN_BOX);
+            DeltaTy_check->deactivate();
+          } // Fl_Check_Button* DeltaTy_check
+          { DeltaTz_check = new Fl_Check_Button(280, 475, 70, 25, "DeltaTz");
+            DeltaTz_check->down_box(FL_DOWN_BOX);
+            DeltaTz_check->deactivate();
+          } // Fl_Check_Button* DeltaTz_check
+          { Hax_check = new Fl_Check_Button(40, 505, 65, 25, "Hax");
+            Hax_check->down_box(FL_DOWN_BOX);
+            Hax_check->deactivate();
+          } // Fl_Check_Button* Hax_check
+          { Hay_check = new Fl_Check_Button(120, 505, 65, 25, "Hay");
+            Hay_check->down_box(FL_DOWN_BOX);
+            Hay_check->deactivate();
+          } // Fl_Check_Button* Hay_check
+          { Za_check = new Fl_Check_Button(200, 505, 65, 25, "Za");
+            Za_check->down_box(FL_DOWN_BOX);
+            Za_check->deactivate();
+          } // Fl_Check_Button* Za_check
+          { mag_para_input = new Fl_Input(220, 440, 240, 28, "Magnetization Parameters :  ");
+            mag_para_input->tooltip("<I0>/<D0>/<I>/<D>");
+            mag_para_input->callback((Fl_Callback*)cb_mag_para_input);
+            mag_para_input->deactivate();
+          } // Fl_Input* mag_para_input
+          mag_group->end();
+        } // Fl_Group* mag_group
+        { mag_data_check = new Fl_Check_Button(40, 405, 170, 30, "Forward Magnetic Data");
+          mag_data_check->down_box(FL_DOWN_BOX);
+          mag_data_check->callback((Fl_Callback*)cb_mag_data_check);
+        } // Fl_Check_Button* mag_data_check
+        { cal_btn = new Fl_Button(40, 545, 420, 28, "Calculate !");
+          cal_btn->callback((Fl_Callback*)cb_cal_btn);
+        } // Fl_Button* cal_btn
+        { noise_para_input = new Fl_Input(260, 370, 200, 28);
+          noise_para_input->tooltip("<noise-mean>/<noise-deviation>");
+          noise_para_input->callback((Fl_Callback*)cb_noise_para_input);
+          noise_para_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+          noise_para_input->deactivate();
+        } // Fl_Input* noise_para_input
+        { noise_check = new Fl_Check_Button(260, 345, 180, 30, "Input Noise Parameters :");
+          noise_check->down_box(FL_DOWN_BOX);
+          noise_check->callback((Fl_Callback*)cb_noise_check);
+        } // Fl_Check_Button* noise_check
+        { res_file_input = new Fl_Input(40, 190, 300, 28, "Input Prefix of Output FIle :");
+          res_file_input->callback((Fl_Callback*)cb_res_file_input);
+          res_file_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Input* res_file_input
+        { res_file_btn = new Fl_Button(360, 190, 100, 28, "Result File");
+          res_file_btn->callback((Fl_Callback*)cb_res_file_btn);
+        } // Fl_Button* res_file_btn
+        { res_out_file_output = new Fl_Output(40, 290, 200, 28, "Prefix of Output File :");
+          res_out_file_output->box(FL_NO_BOX);
+          res_out_file_output->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Output* res_out_file_output
+        forward_tab->end();
+      } // Fl_Group* forward_tab
+      main_tabs->end();
+    } // Fl_Tabs* main_tabs
+    main_window->end();
+  } // Fl_Double_Window* main_window
+  main_window->show(argc, argv);
+  return Fl::run();
+}

fluid_project/gm3d_gui.fl

@@ -0,0 +1,267 @@
+# data file for the Fltk User Interface Designer (fluid)
+version 1.0305
+header_name {.h}
+code_name {.cxx}
+Function {cb_mesh_para_input(Fl_Input*, void*)} {open
+} {}
+
+Function {cb_mesh_file_btn(Fl_Widget*, void*)} {open
+} {}
+
+Function {cb_add_mod_btn(Fl_Widget*, void*)} {open
+} {}
+
+Function {cb_mod_para_file_btn(Fl_Widget*, void*)} {open
+} {}
+
+Function {cb_del_mod_btn(Fl_Widget*, void*)} {open
+} {}
+
+Function {cb_mod_para_brw(Fl_Browser*, void*)} {open
+} {}
+
+Function {cb_mod_ele_input(Fl_Input*, void*)} {open
+} {}
+
+Function {cb_rm_emp_bok_check(Fl_Check_Button*, void*)} {open
+} {}
+
+Function {cb_mod_file_out_btn(Fl_Widget*, void*)} {open
+} {}
+
+Function {cb_build_mod_btn(Fl_Widget*, void*)} {open
+} {}
+
+Function {cb_mod_file_input(Fl_Input*, void*)} {open
+} {}
+
+Function {cb_mod_file_btn(Fl_Widget*, void*)} {open
+} {}
+
+Function {cb_obs_file_input(Fl_Input*, void*)} {open
+} {}
+
+Function {cb_obs_file_btn(Fl_Widget*, void*)} {open
+} {}
+
+Function {cb_res_file_input(Fl_Input*, void*)} {open
+} {}
+
+Function {cb_res_file_btn(Fl_Widget*, void*)} {open
+} {}
+
+Function {cb_nosie_check(Fl_Check_Button*, void*)} {open
+} {}
+
+Function {cb_noise_para_input(Fl_Input*, void*)} {open
+} {}
+
+Function {cb_mag_data_check(Fl_Check_Button*, void*)} {open
+} {}
+
+Function {cb_mag_para_input(Fl_Input*, void*)} {open
+} {}
+
+Function {} {open
+} {
+  Fl_Window main_window {
+    label gm3d open
+    xywh {678 185 500 600} type Double visible
+  } {
+    Fl_Tabs main_tabs {open
+      xywh {10 10 480 580}
+    } {
+      Fl_Group model_tab {
+        label {Build Model} open
+        xywh {10 40 480 550} hide
+      } {
+        Fl_Input mesh_para_input {
+          label {Input Mesh Parameters :}
+          callback cb_mesh_para_input
+          xywh {40 70 300 28} align 5
+        }
+        Fl_Button mesh_file_btn {
+          label {Mesh File}
+          callback cb_mesh_file_btn
+          xywh {360 70 100 28}
+        }
+        Fl_Button mod_para_file_btn {
+          label {Add Model From File}
+          callback cb_mod_para_file_btn
+          xywh {170 160 160 28}
+        }
+        Fl_Input mod_ele_input_build {
+          label {Input Model Element Data Name:}
+          callback cb_mod_ele_input
+          xywh {40 445 215 28} align 5
+        }
+        Fl_Button build_mod_btn {
+          label {Build Model !}
+          callback cb_build_mod_btn
+          xywh {40 545 420 28}
+        }
+        Fl_Output mesh_para_output {
+          label {Mesh Parameters :}
+          xywh {40 120 300 28} box NO_BOX align 5
+        }
+        Fl_Button add_mod_btn {
+          label {Add Model}
+          callback cb_add_mod_btn
+          xywh {40 160 100 28}
+        }
+        Fl_Button del_mod_btn {
+          label {Delete Model}
+          callback cb_del_mod_btn
+          xywh {360 160 100 28} deactivate
+        }
+        Fl_Check_Button rm_emp_bok_check {
+          label {Remove Empty Blocks}
+          callback cb_rm_emp_bok_check
+          xywh {280 447 180 28} down_box DOWN_BOX
+        }
+        Fl_Browser mod_para_brw {
+          callback cb_mod_para_brw
+          xywh {40 205 420 210}
+        }
+        Fl_Button mod_file_out_btn {
+          label {Model File}
+          callback cb_mod_file_out_btn
+          xywh {360 500 100 28}
+        }
+        Fl_Output mod_out_file_output {
+          label {Output File Name :}
+          xywh {40 500 300 28} box NO_BOX align 5
+        }
+      }
+      Fl_Group forward_tab {
+        label {Forward Modeling} open
+        xywh {10 40 480 550}
+      } {
+        Fl_Group grav_group {
+          label {Forward Gravitational Data :} open
+          xywh {40 340 190 58} align 5
+        } {
+          Fl_Check_Button Vz_check {
+            label Vz selected
+            xywh {40 340 50 28} down_box DOWN_BOX
+          }
+          Fl_Check_Button Vzx_check {
+            label Vzx
+            xywh {105 340 50 28} down_box DOWN_BOX
+          }
+          Fl_Check_Button Vzy_check {
+            label Vzy
+            xywh {170 340 50 28} down_box DOWN_BOX
+          }
+          Fl_Check_Button Vzz_check {
+            label Vzz
+            xywh {40 370 50 28} down_box DOWN_BOX
+          }
+        }
+        Fl_Button mod_file_btn {
+          label {Model File}
+          callback cb_mod_file_btn
+          xywh {360 70 100 28}
+        }
+        Fl_Button obs_file_btn {
+          label {Observe File}
+          callback cb_obs_file_btn
+          xywh {360 130 100 28}
+        }
+        Fl_Input mod_file_input {
+          label {Input Model FIle :}
+          callback cb_mod_file_input
+          xywh {40 70 300 28} align 5
+        }
+        Fl_Output mod_file_output {
+          label {Chosen Model File :}
+          xywh {40 240 200 28} box NO_BOX align 5
+        }
+        Fl_Input obs_file_input {
+          label {Input Observe  FIle :}
+          callback cb_obs_file_input
+          xywh {40 130 300 28} align 5
+        }
+        Fl_Output obs_file_output {
+          label {Chosen Observe File :}
+          xywh {260 240 200 28} box NO_BOX align 5
+        }
+        Fl_Input mod_ele_input {
+          label {Model Element Data Name :}
+          callback cb_mod_ele_input
+          xywh {260 315 200 28} align 5
+        }
+        Fl_Group mag_group {open
+          xywh {40 440 420 90}
+        } {
+          Fl_Check_Button DeltaT_check {
+            label DeltaT
+            xywh {40 475 65 25} down_box DOWN_BOX deactivate
+          }
+          Fl_Check_Button DeltaTx_check {
+            label DeltaTx
+            xywh {120 475 70 25} down_box DOWN_BOX deactivate
+          }
+          Fl_Check_Button DeltaTy_check {
+            label DeltaTy
+            xywh {200 475 70 25} down_box DOWN_BOX deactivate
+          }
+          Fl_Check_Button DeltaTz_check {
+            label DeltaTz
+            xywh {280 475 70 25} down_box DOWN_BOX deactivate
+          }
+          Fl_Check_Button Hax_check {
+            label Hax
+            xywh {40 505 65 25} down_box DOWN_BOX deactivate
+          }
+          Fl_Check_Button Hay_check {
+            label Hay
+            xywh {120 505 65 25} down_box DOWN_BOX deactivate
+          }
+          Fl_Check_Button Za_check {
+            label Za
+            xywh {200 505 65 25} down_box DOWN_BOX deactivate
+          }
+          Fl_Input mag_para_input {
+            label {Magnetization Parameters :  }
+            callback cb_mag_para_input
+            tooltip {<I0>/<D0>/<I>/<D>} xywh {220 440 240 28} deactivate
+          }
+        }
+        Fl_Check_Button mag_data_check {
+          label {Forward Magnetic Data}
+          callback cb_mag_data_check
+          xywh {40 405 170 30} down_box DOWN_BOX
+        }
+        Fl_Button cal_btn {
+          label {Calculate !}
+          callback cb_cal_btn
+          xywh {40 545 420 28}
+        }
+        Fl_Input noise_para_input {
+          callback cb_noise_para_input
+          tooltip {<noise-mean>/<nosie-deviation>} xywh {260 370 200 28} align 5 deactivate
+        }
+        Fl_Check_Button noise_check {
+          label {Input Noise Parameters :}
+          callback cb_noise_check
+          xywh {260 345 180 30} down_box DOWN_BOX
+        }
+        Fl_Input res_file_input {
+          label {Input Prefix of Output FIle :}
+          callback cb_res_file_input
+          xywh {40 190 300 28} align 5
+        }
+        Fl_Button res_file_btn {
+          label {Result File}
+          callback cb_res_file_btn
+          xywh {360 190 100 28}
+        }
+        Fl_Output res_out_file_output {
+          label {Prefix of Output File :}
+          xywh {40 290 200 28} box NO_BOX align 5
+        }
+      }
+    }
+  }
+}

fluid_project/gm3d_gui.h

@@ -0,0 +1,82 @@
+// generated by Fast Light User Interface Designer (fluid) version 1.0305
+
+#ifndef _GM3D_GUI_H
+#define _GM3D_GUI_H
+#include <FL/Fl.H>
+#include <FL/Fl_Double_Window.H>
+#include <FL/Fl_Tabs.H>
+#include <FL/Fl_Group.H>
+#include <FL/Fl_Input.H>
+#include <FL/Fl_Button.H>
+#include <FL/Fl_Output.H>
+#include <FL/Fl_Check_Button.H>
+#include <FL/Fl_Browser.H>
+
+extern void cb_mesh_para_input(Fl_Input*, void*);
+extern void cb_mesh_file_btn(Fl_Widget*, void*);
+extern void cb_add_mod_btn(Fl_Widget*, void*);
+extern void cb_mod_para_file_btn(Fl_Widget*, void*);
+extern void cb_del_mod_btn(Fl_Widget*, void*);
+extern void cb_mod_para_brw(Fl_Browser*, void*);
+extern void cb_mod_ele_input(Fl_Input*, void*);
+extern void cb_rm_emp_bok_check(Fl_Check_Button*, void*);
+extern void cb_mod_file_out_btn(Fl_Widget*, void*);
+extern void cb_build_mod_btn(Fl_Widget*, void*);
+extern void cb_mod_file_input(Fl_Input*, void*);
+extern void cb_mod_file_btn(Fl_Widget*, void*);
+extern void cb_obs_file_input(Fl_Input*, void*);
+extern void cb_obs_file_btn(Fl_Widget*, void*);
+extern void cb_res_file_input(Fl_Input*, void*);
+extern void cb_res_file_btn(Fl_Widget*, void*);
+extern void cb_nosie_check(Fl_Check_Button*, void*);
+extern void cb_noise_para_input(Fl_Input*, void*);
+extern void cb_mag_data_check(Fl_Check_Button*, void*);
+extern void cb_mag_para_input(Fl_Input*, void*);
+extern void cb_cal_btn(Fl_Button*, void*);
+extern void cb_noise_check(Fl_Check_Button*, void*);
+
+extern Fl_Double_Window *main_window;
+extern Fl_Tabs *main_tabs;
+extern Fl_Group *model_tab;
+extern Fl_Input *mesh_para_input;
+extern Fl_Button *mesh_file_btn;
+extern Fl_Button *mod_para_file_btn;
+extern Fl_Input *mod_ele_input_build;
+extern Fl_Button *build_mod_btn;
+extern Fl_Output *mesh_para_output;
+extern Fl_Button *add_mod_btn;
+extern Fl_Button *del_mod_btn;
+extern Fl_Check_Button *rm_emp_bok_check;
+extern Fl_Browser *mod_para_brw;
+extern Fl_Button *mod_file_out_btn;
+extern Fl_Output *mod_out_file_output;
+extern Fl_Group *forward_tab;
+extern Fl_Group *grav_group;
+extern Fl_Check_Button *Vz_check;
+extern Fl_Check_Button *Vzx_check;
+extern Fl_Check_Button *Vzy_check;
+extern Fl_Check_Button *Vzz_check;
+extern Fl_Button *mod_file_btn;
+extern Fl_Button *obs_file_btn;
+extern Fl_Input *mod_file_input;
+extern Fl_Output *mod_file_output;
+extern Fl_Input *obs_file_input;
+extern Fl_Output *obs_file_output;
+extern Fl_Input *mod_ele_input;
+extern Fl_Group *mag_group;
+extern Fl_Check_Button *DeltaT_check;
+extern Fl_Check_Button *DeltaTx_check;
+extern Fl_Check_Button *DeltaTy_check;
+extern Fl_Check_Button *DeltaTz_check;
+extern Fl_Check_Button *Hax_check;
+extern Fl_Check_Button *Hay_check;
+extern Fl_Check_Button *Za_check;
+extern Fl_Input *mag_para_input;
+extern Fl_Check_Button *mag_data_check;
+extern Fl_Button *cal_btn;
+extern Fl_Input *noise_para_input;
+extern Fl_Check_Button *noise_check;
+extern Fl_Input *res_file_input;
+extern Fl_Button *res_file_btn;
+extern Fl_Output *res_out_file_output;
+#endif

src/CMakeLists.txt

@@ -0,0 +1,18 @@
+aux_source_directory(. DIR_SRC)
+
+set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin)
+
+find_package(OpenMP)
+if (OPENMP_FOUND)
+	set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+	set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+	set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
+endif()
+# 生成windows和mac的应用程序文件
+add_executable(gm3d_gui WIN32 MACOSX_BUNDLE ${DIR_SRC})
+
+include_directories(/usr/local/include)
+
+find_library(FLTK_LIBRARY fltk HINTS /usr/local/lib)
+
+target_link_libraries(gm3d_gui PUBLIC ${FLTK_LIBRARY})

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_model_gui.cpp

@@ -0,0 +1,182 @@
+#include "gm3d_gui.h"
+
+Fl_Double_Window *add_mod_win=(Fl_Double_Window *)0;
+Fl_Group *mod_type_group=(Fl_Group *)0;
+Fl_Radio_Round_Button *reg_bok_rbtn=(Fl_Radio_Round_Button *)0;
+Fl_Radio_Round_Button *til_bok_rbtn=(Fl_Radio_Round_Button *)0;
+Fl_Radio_Round_Button *sph_rbtn=(Fl_Radio_Round_Button *)0;
+Fl_Radio_Round_Button *int_face_rbtn=(Fl_Radio_Round_Button *)0;
+Fl_Group *val_type_group=(Fl_Group *)0;
+Fl_Radio_Round_Button *app_val_rbtn=(Fl_Radio_Round_Button *)0;
+Fl_Radio_Round_Button *rep_val_rbtn=(Fl_Radio_Round_Button *)0;
+Fl_Radio_Round_Button *era_val_rbtn=(Fl_Radio_Round_Button *)0;
+Fl_Group *agn_part_group=(Fl_Group *)0;
+Fl_Radio_Round_Button *top_val_btn=(Fl_Radio_Round_Button *)0;
+Fl_Radio_Round_Button *bot_val_btn=(Fl_Radio_Round_Button *)0;
+Fl_Input *sig_mod_para_input=(Fl_Input *)0;
+Fl_Return_Button *can_add_btn=(Fl_Return_Button *)0;
+Fl_Button *sig_add_btn=(Fl_Button *)0;
+Fl_Input *mod_val_input=(Fl_Input *)0;
+
+// 声明三个字符串用于保存新建的模型类型和赋值类型
+char add_mod_type[1024];
+char add_val_type[1024];
+char agn_val_type[1024];
+
+void cb_reg_bok_rbtn(Fl_Button*, void*){
+  strcpy(add_mod_type,"regular_block");
+  agn_part_group->deactivate();
+  //printf("%s\n", "Regular block button selected!");
+  return;
+}
+
+void cb_til_bok_rbtn(Fl_Button*, void*){
+  strcpy(add_mod_type,"tilted_block");
+  agn_part_group->deactivate();
+  //printf("%s\n", "Tilted block button selected!");
+  return;
+}
+
+void cb_sph_rbtn(Fl_Button*, void*){
+  strcpy(add_mod_type,"sphere");
+  agn_part_group->deactivate();
+  //printf("%s\n", "Sphere button selected!");
+  return;
+}
+
+void cb_int_face_rbtn(Fl_Button*, void*){
+  strcpy(add_mod_type,"interface");
+  agn_part_group->activate();
+  //printf("%s\n", "Interface button selected!");
+  return;
+}
+
+void cb_app_val_rbtn(Fl_Button*, void*){
+  strcpy(add_val_type,"add");
+  //printf("%s\n", "append button selected!");
+  return;
+}
+
+void cb_rep_val_rbtn(Fl_Button*, void*){
+  strcpy(add_val_type,"replace");
+  //printf("%s\n", "replace button selected!");
+  return;
+}
+
+void cb_era_val_rbtn(Fl_Button*, void*){
+  strcpy(add_val_type,"erase");
+  //printf("%s\n", "Erase button selected!");
+  return;
+}
+
+void cb_top_val_rbtn(Fl_Button*, void*){
+  strcpy(agn_val_type,"top");
+  //printf("%s\n", "Erase button selected!");
+  return;
+}
+
+void cb_bot_val_rbtn(Fl_Button*, void*){
+  strcpy(agn_val_type,"bot");
+  //printf("%s\n", "Erase button selected!");
+  return;
+}
+
+// 声明一个字符串用于保存新建的模型参数并添加至列表中
+// 模型参数的排列为<模型类型> <赋值类型> <物理参数字符串> <几何参数字符串>
+char add_mod_para[1024];
+
+void cb_sig_add_btn(Fl_Widget*, void*){
+  strcpy(add_mod_para,add_mod_type);
+  strcat(add_mod_para," ");
+  strcat(add_mod_para,add_val_type);
+  if (!strcmp(add_mod_type,"interface"))
+  {
+    strcat(add_mod_para,"/");
+    strcat(add_mod_para,agn_val_type);
+  }
+  strcat(add_mod_para," ");
+  strcat(add_mod_para,mod_val_input->value());
+  strcat(add_mod_para," ");
+  strcat(add_mod_para,sig_mod_para_input->value());
+  mod_para_brw->add(add_mod_para);
+  return;
+}
+//取消添加模型
+void cb_can_add_btn(Fl_Widget*, void*){
+	add_mod_win->hide();
+	return;
+}
+
+void cb_add_mod_btn(Fl_Widget*, void*)
+{
+  { add_mod_win = new Fl_Double_Window(315, 300, "Add model (gm3d)");
+    { mod_type_group = new Fl_Group(20, 25, 250, 53, "Model Type :");
+      mod_type_group->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+      { reg_bok_rbtn = new Fl_Radio_Round_Button(20, 25, 110, 28, "Regular Block");
+        reg_bok_rbtn->down_box(FL_ROUND_DOWN_BOX);
+        reg_bok_rbtn->callback((Fl_Callback*)cb_reg_bok_rbtn);
+        reg_bok_rbtn->tooltip("model parameter:\n<xmin>/<xmax>/<ymin>/<ymax>/<zmin>/<zmax>");
+      } // Fl_Radio_Round_Button* reg_bok_rbtn
+      { til_bok_rbtn = new Fl_Radio_Round_Button(160, 25, 100, 28, "Tilted Block");
+        til_bok_rbtn->down_box(FL_ROUND_DOWN_BOX);
+        til_bok_rbtn->callback((Fl_Callback*)cb_til_bok_rbtn);
+        til_bok_rbtn->tooltip("model parameter:\n<xmin_z>/<xmax_z>/<ymin_z>/<ymax_z>/<zmin>/<xmin_Z>/<xmax_Z>/<ymin_Z>/<ymax_Z>/<zmax>");
+      } // Fl_Radio_Round_Button* til_bok_rbtn
+      { sph_rbtn = new Fl_Radio_Round_Button(20, 50, 70, 28, "Sphere");
+        sph_rbtn->down_box(FL_ROUND_DOWN_BOX);
+        sph_rbtn->callback((Fl_Callback*)cb_sph_rbtn);
+        sph_rbtn->tooltip("model parameter:\n<x_c>/<y_c>/<z_c>/<x_radius>/<y_radius>/<z_radius>");
+      } // Fl_Radio_Round_Button* sph_rbtn
+      { int_face_rbtn = new Fl_Radio_Round_Button(160, 50, 80, 28, "Interface");
+        int_face_rbtn->down_box(FL_ROUND_DOWN_BOX);
+        int_face_rbtn->callback((Fl_Callback*)cb_int_face_rbtn);
+        int_face_rbtn->tooltip("model parameter:\n<grid-filename>");
+      } // Fl_Radio_Round_Button* int_face_rbtn
+      mod_type_group->end();
+    } // Fl_Group* mod_type_group
+    { val_type_group = new Fl_Group(20, 100, 250, 28, "Value Type :");
+      val_type_group->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+      { app_val_rbtn = new Fl_Radio_Round_Button(20, 100, 80, 28, "Append");
+        app_val_rbtn->down_box(FL_ROUND_DOWN_BOX);
+        app_val_rbtn->callback((Fl_Callback*)cb_app_val_rbtn);
+      } // Fl_Radio_Round_Button* app_val_rbtn
+      { rep_val_rbtn = new Fl_Radio_Round_Button(110, 100, 80, 28, "Replace");
+        rep_val_rbtn->down_box(FL_ROUND_DOWN_BOX);
+        rep_val_rbtn->callback((Fl_Callback*)cb_rep_val_rbtn);
+      } // Fl_Radio_Round_Button* rep_val_rbtn
+      { era_val_rbtn = new Fl_Radio_Round_Button(200, 100, 60, 28, "Erase");
+        era_val_rbtn->down_box(FL_ROUND_DOWN_BOX);
+        era_val_rbtn->callback((Fl_Callback*)cb_era_val_rbtn);
+      } // Fl_Radio_Round_Button* era_val_rbtn
+      val_type_group->end();
+    } // Fl_Group* val_type_group
+    { agn_part_group = new Fl_Group(20, 150, 160, 28, "Assign Part :");
+      agn_part_group->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+      { top_val_btn = new Fl_Radio_Round_Button(20, 150, 50, 28, "Top");
+        top_val_btn->down_box(FL_ROUND_DOWN_BOX);
+        top_val_btn->callback((Fl_Callback*)cb_top_val_rbtn);
+      } // Fl_Radio_Round_Button* top_val_btn
+      { bot_val_btn = new Fl_Radio_Round_Button(110, 150, 70, 28, "Bottom");
+        bot_val_btn->down_box(FL_ROUND_DOWN_BOX);
+        bot_val_btn->callback((Fl_Callback*)cb_bot_val_rbtn);
+      } // Fl_Radio_Round_Button* bot_val_btn
+      agn_part_group->deactivate();
+      agn_part_group->end();
+    } // Fl_Group* agn_part_group
+    { sig_mod_para_input = new Fl_Input(20, 200, 275, 28, "Model Parameter :");
+      sig_mod_para_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+    } // Fl_Input* sig_mod_para_input
+    { can_add_btn = new Fl_Return_Button(210, 255, 85, 28, "Cancel");
+      can_add_btn->callback((Fl_Callback*)cb_can_add_btn);
+    } // Fl_Return_Button* can_add_btn
+    { sig_add_btn = new Fl_Button(140, 255, 60, 28, "Add");
+      sig_add_btn->callback((Fl_Callback*)cb_sig_add_btn);
+    } // Fl_Button* sig_add_btn
+    { mod_val_input = new Fl_Input(20, 255, 110, 28, "Model Value :");
+      mod_val_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+    } // Fl_Input* mod_val_input
+    add_mod_win->end();
+    add_mod_win->show();
+  } // Fl_Double_Window* add_mod_win
+  return;
+}

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_models_gui.cpp

@@ -0,0 +1,23 @@
+#include "gm3d.h"
+
+void GM3D::AddModels_GUI(){
+	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;
+}

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,88 @@
+#include "gm3d.h"
+
+int GM3D::BuildRegularGrid(char* space_para){
+	cpoint temp_cp;
+	cube temp_cu;
+	string temp_str;
+	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;
+
+	ifstream infile;
+	infile.open(space_para);
+	if (!infile) {
+		if (9 != sscanf(space_para,"%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf",
+					&xs,&dx,&xe,&ys,&dy,&ye,&zs,&dz,&ze)){
+			return -1;
+		}
+	}
+	else {
+		getline(infile,temp_str);
+		if (3 != sscanf(temp_str.c_str(),"%lf %lf %lf",&xs,&dx,&xe)) return -1;
+
+		getline(infile,temp_str);
+		if (3 != sscanf(temp_str.c_str(),"%lf %lf %lf",&ys,&dy,&ye)) return -1;
+
+		getline(infile,temp_str);
+		if (3 != sscanf(temp_str.c_str(),"%lf %lf %lf",&zs,&dz,&ze)) return -1;
+
+		infile.close();
+	}
+
+	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()){
+		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/forward_delta_t.cpp

@@ -0,0 +1,136 @@
+#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;
+
+	//初始化正演值和不确定度
+	for (int i = 0; i < obs_num_; i++){
+		obs_p_[i].val = obs_p_[i].dev = 0.0;
+	}
+
+	//确定噪声水平
+	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,113 @@
+#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;
+
+	//初始化正演值和不确定度
+	for (int i = 0; i < obs_num_; i++){
+		obs_p_[i].val = obs_p_[i].dev = 0.0;
+	}
+
+	//确定噪声水平
+	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,113 @@
+#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;
+
+	//初始化正演值和不确定度
+	for (int i = 0; i < obs_num_; i++){
+		obs_p_[i].val = obs_p_[i].dev = 0.0;
+	}
+
+	//确定噪声水平
+	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,124 @@
+#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;
+
+	//初始化正演值和不确定度
+	for (int i = 0; i < obs_num_; i++){
+		obs_p_[i].val = obs_p_[i].dev = 0.0;
+	}
+
+	//确定噪声水平
+	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,108 @@
+#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;
+
+	//初始化正演值和不确定度
+	for (int i = 0; i < obs_num_; i++){
+		obs_p_[i].val = obs_p_[i].dev = 0.0;
+	}
+
+	//确定噪声水平
+	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,108 @@
+#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;
+
+	//初始化正演值和不确定度
+	for (int i = 0; i < obs_num_; i++){
+		obs_p_[i].val = obs_p_[i].dev = 0.0;
+	}
+
+	//确定噪声水平
+	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,69 @@
+#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;
+
+	//初始化正演值和不确定度
+	for (int i = 0; i < obs_num_; i++){
+		obs_p_[i].val = obs_p_[i].dev = 0.0;
+	}
+
+	//确定噪声水平
+	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,68 @@
+#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;
+
+	//初始化正演值和不确定度
+	for (int i = 0; i < obs_num_; i++){
+		obs_p_[i].val = obs_p_[i].dev = 0.0;
+	}
+
+	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,68 @@
+#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;
+
+	//初始化正演值和不确定度
+	for (int i = 0; i < obs_num_; i++){
+		obs_p_[i].val = obs_p_[i].dev = 0.0;
+	}
+
+	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,77 @@
+#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;
+
+	//初始化正演值和不确定度
+	for (int i = 0; i < obs_num_; i++){
+		obs_p_[i].val = obs_p_[i].dev = 0.0;
+	}
+
+	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,108 @@
+#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;
+
+	//初始化正演值和不确定度
+	for (int i = 0; i < obs_num_; i++){
+		obs_p_[i].val = obs_p_[i].dev = 0.0;
+	}
+
+	//确定噪声水平
+	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/get_model_list.cpp

@@ -0,0 +1,6 @@
+#include "gm3d.h"
+
+void GM3D::get_model_list(modelistArray input_list){
+	model_list_ = input_list;
+	return;
+}

src/gm3d.h

@@ -0,0 +1,61 @@
+#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*); //读取模型块体参数文件
+	void AddModels_GUI(); //添加模型块体参数
+	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*);
+	//gm3d_gui函数
+	void get_model_list(modelistArray);
+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/gm3d_gui.cpp

@@ -0,0 +1,539 @@
+#include "gm3d_gui.h"
+#include "gm3d.h"
+
+Fl_Double_Window *main_window=(Fl_Double_Window *)0;
+Fl_Tabs *main_tabs=(Fl_Tabs *)0;
+Fl_Group *model_tab=(Fl_Group *)0;
+Fl_Input *mesh_para_input=(Fl_Input *)0;
+Fl_Button *mesh_file_btn=(Fl_Button *)0;
+Fl_Button *mod_para_file_btn=(Fl_Button *)0;
+Fl_Input *mod_ele_input_build=(Fl_Input *)0;
+Fl_Button *build_mod_btn=(Fl_Button *)0;
+Fl_Output *mesh_para_output=(Fl_Output *)0;
+Fl_Button *add_mod_btn=(Fl_Button *)0;
+Fl_Button *del_mod_btn=(Fl_Button *)0;
+Fl_Check_Button *rm_emp_bok_check=(Fl_Check_Button *)0;
+Fl_Browser *mod_para_brw=(Fl_Browser *)0;
+Fl_Button *mod_file_out_btn=(Fl_Button *)0;
+Fl_Output *mod_out_file_output=(Fl_Output *)0;
+Fl_Group *forward_tab=(Fl_Group *)0;
+Fl_Group *grav_group=(Fl_Group *)0;
+Fl_Check_Button *Vz_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Vzx_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Vzy_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Vzz_check=(Fl_Check_Button *)0;
+Fl_Button *mod_file_btn=(Fl_Button *)0;
+Fl_Button *obs_file_btn=(Fl_Button *)0;
+Fl_Input *mod_file_input=(Fl_Input *)0;
+Fl_Output *mod_file_output=(Fl_Output *)0;
+Fl_Input *obs_file_input=(Fl_Input *)0;
+Fl_Output *obs_file_output=(Fl_Output *)0;
+Fl_Input *mod_ele_input=(Fl_Input *)0;
+Fl_Group *mag_group=(Fl_Group *)0;
+Fl_Check_Button *DeltaT_check=(Fl_Check_Button *)0;
+Fl_Check_Button *DeltaTx_check=(Fl_Check_Button *)0;
+Fl_Check_Button *DeltaTy_check=(Fl_Check_Button *)0;
+Fl_Check_Button *DeltaTz_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Hax_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Hay_check=(Fl_Check_Button *)0;
+Fl_Check_Button *Za_check=(Fl_Check_Button *)0;
+Fl_Input *mag_para_input=(Fl_Input *)0;
+Fl_Check_Button *mag_data_check=(Fl_Check_Button *)0;
+Fl_Button *cal_btn=(Fl_Button *)0;
+Fl_Input *noise_para_input=(Fl_Input *)0;
+Fl_Check_Button *noise_check=(Fl_Check_Button *)0;
+Fl_Input *res_file_input=(Fl_Input *)0;
+Fl_Button *res_file_btn=(Fl_Button *)0;
+Fl_Output *res_file_output=(Fl_Output *)0;
+
+char mesh_filename[1024];
+char out_msh_filename[1024];
+
+void cb_mesh_para_input(Fl_Input*, void*){
+	strcpy(mesh_filename,mesh_para_input->value());
+	mesh_para_output->value(mesh_para_input->value());
+	return;
+}
+
+void cb_mesh_file_btn(Fl_Widget*, void*){
+  // Create the file chooser, and show it
+  // directory, filter, chooser type, title
+  Fl_File_Chooser chooser(".", "*.txt, *.dat, *.xyz", Fl_File_Chooser::SINGLE, "Choose a mesh file.");
+  chooser.show();
+
+  // Block until user picks something.
+  while(chooser.shown()) { Fl::wait(); }
+
+  //从文件的绝对路径中分离文件名并从输出框输出
+  if (chooser.value() != NULL)
+  {
+    strcpy(mesh_filename,chooser.value());
+
+    std::string abs_filename = chooser.value();
+    int pos = abs_filename.find_last_of("/");
+
+    std::stringstream temp_ss;
+    temp_ss << abs_filename.substr(pos+1);
+    const char* disp_filename = temp_ss.str().c_str();
+    mesh_para_output->value(disp_filename);
+  }
+  else
+  {
+    mesh_para_output->value("Unset");
+  }
+  return;
+}
+
+void cb_mod_para_file_btn(Fl_Widget*, void*){
+  // Create the file chooser, and show it
+  Fl_File_Chooser chooser(".", "*.txt, *.dat, *.xyz", Fl_File_Chooser::SINGLE, "Choose a model list file.");
+  chooser.show();
+
+  // Block until user picks something.
+  while(chooser.shown()) { Fl::wait(); }
+
+  //下面我们读入文件并添加模型参数到列表显示
+  std::ifstream modpara_in;
+  modpara_in.open(chooser.value());
+  if (!modpara_in)
+  {
+    fl_message("Error! Fail to open file: %s",chooser.value());
+    return;
+  }
+  else{
+    std::string temp_str;
+    std::stringstream temp_ss;
+
+    while(getline(modpara_in,temp_str)){
+      if (*(temp_str.begin()) == '#') continue;
+      else{
+        temp_ss.clear(); temp_ss.str("");
+        temp_ss << temp_str;
+        mod_para_brw->add(temp_ss.str().c_str());
+      }
+    }
+    modpara_in.close();
+  }
+  return;
+}
+
+void cb_del_mod_btn(Fl_Widget*, void*){
+  int index = mod_para_brw->value();
+  mod_para_brw->remove(index);
+  del_mod_btn->deactivate();
+}
+
+void cb_mod_para_brw(Fl_Browser*, void*){
+  for ( int t=1; t<=mod_para_brw->size(); t++ ) {
+    if ( mod_para_brw->selected(t) ) {
+      //printf("%d) '%s'\n", t, mod_para_brw->text(t));
+      del_mod_btn->activate();
+    }
+  }
+  return;
+}
+
+bool remove_null = false;
+void cb_rm_emp_bok_check(Fl_Check_Button*, void*){
+  if (rm_emp_bok_check->value() == 1) remove_null = true;
+  else if (rm_emp_bok_check->value() == 0) remove_null = false;
+  return;
+}
+
+void cb_mod_file_out_btn(Fl_Widget*, void*){
+  // Create the file chooser, and show it
+  // directory, filter, chooser type, title
+  Fl_File_Chooser chooser(".", "*.msh", Fl_File_Chooser::CREATE, "Create a new model file.");
+  chooser.show();
+
+  // Block until user picks something.
+  while(chooser.shown()) { Fl::wait(); }
+
+  //从文件的绝对路径中分离文件名并从输出框输出
+  if (chooser.value() != NULL)
+  {
+    strcpy(out_msh_filename,chooser.value());
+
+    std::string abs_filename = chooser.value();
+    int pos = abs_filename.find_last_of("/");
+
+    std::stringstream temp_ss;
+    temp_ss << abs_filename.substr(pos+1);
+    const char* disp_filename = temp_ss.str().c_str();
+    mod_out_file_output->value(disp_filename);
+  }
+  else
+  {
+    mod_out_file_output->value("Unset");
+  }
+  return;
+}
+
+void cb_build_mod_btn(Fl_Widget*, void*){
+  //首先构建一个GM3D实例
+  GM3D gm3d_instance;
+  //设置初始化变量
+  char dimension[1024] = "10/20/990/10/20/990/10/20/490";
+  char out_mshname[1024] = "Untitled.msh";
+  char elename[1024] = "Untitled";
+
+  //将mesh_para_output拷贝到dimension
+  strcpy(dimension, mesh_filename);
+  //std::cout << dimension << std::endl;
+  //将mod_para_brw中的模型参数列表拷贝到一个临时参数列表
+  std::stringstream temp_ss;
+  modelist temp_list;
+  modelistArray brw_model_list;
+
+  for ( int t=1; t<=mod_para_brw->size(); t++ ) {
+      temp_ss.clear(); temp_ss.str("");
+      temp_ss << mod_para_brw->text(t);
+      if(temp_ss >> temp_list.mod_type >> temp_list.val_type >> temp_list.mod_value >> temp_list.mod_para){
+        brw_model_list.push_back(temp_list);
+      }
+      else{
+        fl_message("Wrong model parameter: %s. Skipped !",mod_para_brw->text(t));
+      }
+  }
+
+  //将brw_model_list拷贝到GM3D中
+  gm3d_instance.get_model_list(brw_model_list);
+  //将mod_out_file_input拷贝到out_mshname
+  strcpy(out_mshname, out_msh_filename);
+  if (!strcmp(out_mshname,"")){
+    fl_message("Output file's name can't be empty !");
+    return;
+  }
+  //将mod_ele_input_build拷贝到elename
+  strcpy(elename, mod_ele_input_build->value());
+  if (!strcmp(elename,"")){
+    fl_message("Element data's name can't be empty !");
+    return;
+  }
+  //构建模型网络
+  if(gm3d_instance.BuildRegularGrid(dimension)){
+    fl_message("Mesh Parameters Load Error !");
+    return;
+  }
+  gm3d_instance.AddModels_GUI();
+  gm3d_instance.RegisteredOuput(remove_null);
+  if (gm3d_instance.OutMshFile(out_mshname,elename)) return;
+
+  //这里需要一个输出信息窗口
+  fl_message("Model Construction Completed !");
+  return;
+}
+
+/***********************************以下是正演计算GUI函数******************************/
+char in_msh_filename[1024];
+char in_obs_filename[1024];
+char out_res_filename[1024];
+
+void cb_mod_file_input(Fl_Input*, void*){
+	strcpy(in_msh_filename,mod_file_input->value());
+	mod_file_output->value(mod_file_input->value());
+	return;
+}
+
+void cb_mod_file_btn(Fl_Widget*, void*){
+  // Create the file chooser, and show it
+  // directory, filter, chooser type, title
+  Fl_File_Chooser chooser(".", "*.msh", Fl_File_Chooser::SINGLE, "Create a model file.");
+  chooser.show();
+
+  // Block until user picks something.
+  while(chooser.shown()) { Fl::wait(); }
+
+  //从文件的绝对路径中分离文件名并从输出框输出
+  if (chooser.value() != NULL)
+  {
+    strcpy(in_msh_filename,chooser.value());
+
+    std::string abs_filename = chooser.value();
+    int pos = abs_filename.find_last_of("/");
+
+    std::stringstream temp_ss;
+    temp_ss << abs_filename.substr(pos+1);
+    const char* disp_filename = temp_ss.str().c_str();
+    mod_file_output->value(disp_filename);
+  }
+  else
+  {
+    mod_file_output->value("Unset");
+  }
+  return;
+}
+
+void cb_obs_file_input(Fl_Input*, void*){
+	strcpy(in_obs_filename,obs_file_input->value());
+	obs_file_output->value(obs_file_input->value());
+	return;
+}
+
+void cb_obs_file_btn(Fl_Widget*, void*){
+  // Create the file chooser, and show it
+  // directory, filter, chooser type, title
+  Fl_File_Chooser chooser(".", "*.txt, *.dat, *.xyz", Fl_File_Chooser::SINGLE, "Create observe file.");
+  chooser.show();
+
+  // Block until user picks something.
+  while(chooser.shown()) { Fl::wait(); }
+
+  //从文件的绝对路径中分离文件名并从输出框输出
+  if (chooser.value() != NULL)
+  {
+    strcpy(in_obs_filename,chooser.value());
+
+    std::string abs_filename = chooser.value();
+    int pos = abs_filename.find_last_of("/");
+
+    std::stringstream temp_ss;
+    temp_ss << abs_filename.substr(pos+1);
+    const char* disp_filename = temp_ss.str().c_str();
+    obs_file_output->value(disp_filename);
+  }
+  else
+  {
+    obs_file_output->value("Unset");
+  }
+  return;
+}
+
+void cb_res_file_input(Fl_Input*, void*){
+	strcpy(out_res_filename,res_file_input->value());
+	res_file_output->value(res_file_input->value());
+	return;
+}
+
+void cb_res_file_btn(Fl_Widget*, void*){
+  // Create the file chooser, and show it
+  // directory, filter, chooser type, title
+  Fl_File_Chooser chooser(".", "*.txt, *.dat", Fl_File_Chooser::CREATE, "Create a new observe file.");
+  chooser.show();
+
+  // Block until user picks something.
+  while(chooser.shown()) { Fl::wait(); }
+
+  //从文件的绝对路径中分离文件名并从输出框输出
+  if (chooser.value() != NULL)
+  {
+    strcpy(out_res_filename,chooser.value());
+
+    std::string abs_filename = chooser.value();
+    int pos = abs_filename.find_last_of("/");
+
+    std::stringstream temp_ss;
+    temp_ss << abs_filename.substr(pos+1);
+    const char* disp_filename = temp_ss.str().c_str();
+    res_file_output->value(disp_filename);
+  }
+  else
+  {
+    res_file_output->value("Unset");
+  }
+  return;
+}
+
+void cb_noise_check(Fl_Check_Button*, void*){
+  if (noise_check->value() == 1){
+    noise_para_input->activate();
+  }
+  else if (noise_check->value() == 0){
+    noise_para_input->deactivate();
+  }
+  return;
+}
+
+void cb_mag_data_check(Fl_Check_Button*, void*){
+  if (mag_data_check->value() == 1){
+    mag_para_input->activate();
+    DeltaT_check->activate();
+    DeltaTx_check->activate();
+    DeltaTy_check->activate();
+    DeltaTz_check->activate();
+    Hax_check->activate();
+    Hay_check->activate();
+    Za_check->activate();
+  }
+  else if (mag_data_check->value() == 0){
+    mag_para_input->deactivate();
+    DeltaT_check->deactivate();
+    DeltaTx_check->deactivate();
+    DeltaTy_check->deactivate();
+    DeltaTz_check->deactivate();
+    Hax_check->deactivate();
+    Hay_check->deactivate();
+    Za_check->deactivate();
+  }
+  return;
+}
+
+void cb_cal_btn(Fl_Button*, void*){
+  //首先构建一个GM3D实例
+  GM3D gm3d_instance;
+  //设置初始化变量
+  char in_mshname[1024] = "Untitled.msh";
+  char in_obspara[1024] = "Untitled.txt";
+  char res_outfile[1024] = "Untitled";
+  char res_outfile_full[1024];
+  char ele_name[1024] = "Untitled";
+  char noise_para[1024] = "0.0/0.0";
+  char mag_para[1024] = "90.0/0.0/90.0/0.0";
+
+  strcpy(in_mshname,in_msh_filename);
+  strcpy(in_obspara,in_obs_filename);
+  strcpy(res_outfile,out_res_filename);
+  strcpy(ele_name,mod_ele_input->value());
+  strcpy(noise_para,noise_para_input->value());
+  strcpy(mag_para,mag_para_input->value());
+
+  if (gm3d_instance.ReadModel(in_mshname,ele_name))
+  {
+    fl_message("Model Load Error !");
+    return;
+  }
+
+  if (gm3d_instance.InitObs(in_obspara))
+  {
+    fl_message("Observe File Load Error !");
+    return;
+  }
+
+  if (Vz_check->value() == 1)
+  {
+    gm3d_instance.ForwardVz(noise_para);
+    strcpy(res_outfile_full,res_outfile);
+    strcat(res_outfile_full,"_Vz.txt");
+    if (gm3d_instance.OutObs(res_outfile_full))
+    {
+      fl_message("Observe File Output Error !");
+      return;
+    }
+  }
+
+  if (Vzx_check->value() == 1)
+  {
+    gm3d_instance.ForwardVzx(noise_para);
+    strcpy(res_outfile_full,res_outfile);
+    strcat(res_outfile_full,"_Vzx.txt");
+    if (gm3d_instance.OutObs(res_outfile_full))
+    {
+      fl_message("Observe File Output Error !");
+      return;
+    }
+  }
+
+  if (Vzy_check->value() == 1)
+  {
+    gm3d_instance.ForwardVzy(noise_para);
+    strcpy(res_outfile_full,res_outfile);
+    strcat(res_outfile_full,"_Vzy.txt");
+    if (gm3d_instance.OutObs(res_outfile_full))
+    {
+      fl_message("Observe File Output Error !");
+      return;
+    }
+  }
+
+  if (Vzz_check->value() == 1)
+  {
+    gm3d_instance.ForwardVzz(noise_para);
+    strcpy(res_outfile_full,res_outfile);
+    strcat(res_outfile_full,"_Vzz.txt");
+    if (gm3d_instance.OutObs(res_outfile_full))
+    {
+      fl_message("Observe File Output Error !");
+      return;
+    }
+  }
+
+  if (mag_data_check->value() == 1)
+  {
+    if (DeltaT_check->value() == 1)
+    {
+      gm3d_instance.ForwardDeltaT(noise_para,mag_para);
+      strcpy(res_outfile_full,res_outfile);
+      strcat(res_outfile_full,"_DeltaT.txt");
+      if (gm3d_instance.OutObs(res_outfile_full))
+      {
+        fl_message("Observe File Output Error !");
+        return;
+      }
+    }
+
+    if (DeltaTx_check->value() == 1)
+    {
+      gm3d_instance.ForwardDeltaTx(noise_para,mag_para);
+      strcpy(res_outfile_full,res_outfile);
+      strcat(res_outfile_full,"_DeltaTx.txt");
+      if (gm3d_instance.OutObs(res_outfile_full))
+      {
+        fl_message("Observe File Output Error !");
+        return;
+      }
+    }
+
+    if (DeltaTy_check->value() == 1)
+    {
+      gm3d_instance.ForwardDeltaTy(noise_para,mag_para);
+      strcpy(res_outfile_full,res_outfile);
+      strcat(res_outfile_full,"_DeltaTy.txt");
+      if (gm3d_instance.OutObs(res_outfile_full))
+      {
+        fl_message("Observe File Output Error !");
+        return;
+      }
+    }
+
+    if (DeltaTz_check->value() == 1)
+    {
+      gm3d_instance.ForwardDeltaTz(noise_para,mag_para);
+      strcpy(res_outfile_full,res_outfile);
+      strcat(res_outfile_full,"_DeltaTz.txt");
+      if (gm3d_instance.OutObs(res_outfile_full))
+      {
+        fl_message("Observe File Output Error !");
+        return;
+      }
+    }
+
+    if (Hax_check->value() == 1)
+    {
+      gm3d_instance.ForwardHax(noise_para,mag_para);
+      strcpy(res_outfile_full,res_outfile);
+      strcat(res_outfile_full,"_Hax.txt");
+      if (gm3d_instance.OutObs(res_outfile_full))
+      {
+        fl_message("Observe File Output Error !");
+        return;
+      }
+    }
+
+    if (Hay_check->value() == 1)
+    {
+      gm3d_instance.ForwardHay(noise_para,mag_para);
+      strcpy(res_outfile_full,res_outfile);
+      strcat(res_outfile_full,"_Hay.txt");
+      if (gm3d_instance.OutObs(res_outfile_full))
+      {
+        fl_message("Observe File Output Error !");
+        return;
+      }
+    }
+
+    if (Za_check->value() == 1)
+    {
+      gm3d_instance.ForwardZa(noise_para,mag_para);
+      strcpy(res_outfile_full,res_outfile);
+      strcat(res_outfile_full,"_Za.txt");
+      if (gm3d_instance.OutObs(res_outfile_full))
+      {
+        fl_message("Observe File Output Error !");
+        return;
+      }
+    }
+  }
+
+  fl_message("Forward calculation completed !");
+  return;
+}

src/gm3d_gui.h

@@ -0,0 +1,115 @@
+// generated by Fast Light User Interface Designer (fluid) version 1.0305
+
+#ifndef _GM3D_GUI_H
+#define _GM3D_GUI_H
+#include <FL/Fl.H>
+#include <FL/Fl_Double_Window.H>
+#include <FL/Fl_Tabs.H>
+#include <FL/Fl_Group.H>
+#include <FL/Fl_Input.H>
+#include <FL/Fl_Button.H>
+#include <FL/Fl_Output.H>
+#include <FL/Fl_Check_Button.H>
+#include <FL/Fl_Radio_Round_Button.H>
+#include <FL/Fl_Browser.H>
+#include <FL/Fl_File_Chooser.H>
+
+#include <iostream>
+#include <string.h>
+#include <sstream>
+#include <fstream>
+
+extern void cb_mesh_para_input(Fl_Input*, void*);
+extern void cb_mesh_file_btn(Fl_Widget*, void*);
+extern void cb_add_mod_btn(Fl_Widget*, void*);
+extern void cb_mod_para_file_btn(Fl_Widget*, void*);
+extern void cb_del_mod_btn(Fl_Widget*, void*);
+extern void cb_mod_para_brw(Fl_Browser*, void*);
+extern void cb_rm_emp_bok_check(Fl_Check_Button*, void*);
+extern void cb_mod_file_out_btn(Fl_Widget*, void*);
+extern void cb_build_mod_btn(Fl_Widget*, void*);
+extern void cb_mod_file_input(Fl_Input*, void*);
+extern void cb_mod_file_btn(Fl_Widget*, void*);
+extern void cb_obs_file_input(Fl_Input*, void*);
+extern void cb_obs_file_btn(Fl_Widget*, void*);
+extern void cb_res_file_input(Fl_Input*, void*);
+extern void cb_res_file_btn(Fl_Widget*, void*);
+extern void cb_noise_check(Fl_Check_Button*, void*);
+extern void cb_mag_data_check(Fl_Check_Button*, void*);
+extern void cb_cal_btn(Fl_Button*, void*);
+
+extern void cb_reg_bok_rbtn(Fl_Button*, void*);
+extern void cb_til_bok_rbtn(Fl_Button*, void*);
+extern void cb_sph_rbtn(Fl_Button*, void*);
+extern void cb_int_face_rbtn(Fl_Button*, void*);
+extern void cb_app_val_rbtn(Fl_Button*, void*);
+extern void cb_rep_val_rbtn(Fl_Button*, void*);
+extern void cb_era_val_rbtn(Fl_Button*, void*);
+extern void cb_top_val_rbtn(Fl_Button*, void*);
+extern void cb_bot_val_rbtn(Fl_Button*, void*);
+extern void cb_sig_add_btn(Fl_Widget*, void*);
+extern void cb_can_add_btn(Fl_Widget*, void*);
+
+extern Fl_Double_Window *main_window;
+extern Fl_Tabs *main_tabs;
+extern Fl_Group *model_tab;
+extern Fl_Input *mesh_para_input;
+extern Fl_Button *mesh_file_btn;
+extern Fl_Button *mod_para_file_btn;
+extern Fl_Input *mod_ele_input_build;
+extern Fl_Button *build_mod_btn;
+extern Fl_Output *mesh_para_output;
+extern Fl_Button *add_mod_btn;
+extern Fl_Button *del_mod_btn;
+extern Fl_Check_Button *rm_emp_bok_check;
+extern Fl_Browser *mod_para_brw;
+extern Fl_Button *mod_file_out_btn;
+extern Fl_Output *mod_out_file_output;
+extern Fl_Group *forward_tab;
+extern Fl_Group *grav_group;
+extern Fl_Check_Button *Vz_check;
+extern Fl_Check_Button *Vzx_check;
+extern Fl_Check_Button *Vzy_check;
+extern Fl_Check_Button *Vzz_check;
+extern Fl_Button *mod_file_btn;
+extern Fl_Button *obs_file_btn;
+extern Fl_Input *mod_file_input;
+extern Fl_Output *mod_file_output;
+extern Fl_Input *obs_file_input;
+extern Fl_Output *obs_file_output;
+extern Fl_Input *mod_ele_input;
+extern Fl_Group *mag_group;
+extern Fl_Check_Button *DeltaT_check;
+extern Fl_Check_Button *DeltaTx_check;
+extern Fl_Check_Button *DeltaTy_check;
+extern Fl_Check_Button *DeltaTz_check;
+extern Fl_Check_Button *Hax_check;
+extern Fl_Check_Button *Hay_check;
+extern Fl_Check_Button *Za_check;
+extern Fl_Input *mag_para_input;
+extern Fl_Check_Button *mag_data_check;
+extern Fl_Button *cal_btn;
+extern Fl_Input *noise_para_input;
+extern Fl_Check_Button *noise_check;
+extern Fl_Input *res_file_input;
+extern Fl_Button *res_file_btn;
+extern Fl_Output *res_file_output;
+
+extern Fl_Double_Window *add_mod_win;
+extern Fl_Group *mod_type_group;
+extern Fl_Radio_Round_Button *reg_bok_rbtn;
+extern Fl_Radio_Round_Button *til_bok_rbtn;
+extern Fl_Radio_Round_Button *sph_rbtn;
+extern Fl_Radio_Round_Button *int_face_rbtn;
+extern Fl_Group *val_type_group;
+extern Fl_Radio_Round_Button *app_val_rbtn;
+extern Fl_Radio_Round_Button *rep_val_rbtn;
+extern Fl_Radio_Round_Button *era_val_rbtn;
+extern Fl_Group *agn_part_group;
+extern Fl_Radio_Round_Button *top_val_btn;
+extern Fl_Radio_Round_Button *bot_val_btn;
+extern Fl_Input *sig_mod_para_input;
+extern Fl_Return_Button *can_add_btn;
+extern Fl_Button *sig_add_btn;
+extern Fl_Input *mod_val_input;
+#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,177 @@
+// generated by Fast Light User Interface Designer (fluid) version 1.0305
+
+#include "gm3d_gui.h"
+
+int main(int argc, char **argv) {
+  { main_window = new Fl_Double_Window(500, 600, "gm3d");
+    { main_tabs = new Fl_Tabs(10, 10, 480, 580);
+      { model_tab = new Fl_Group(10, 40, 480, 550, "Build Model");
+        //model_tab->hide();
+        { mesh_para_input = new Fl_Input(40, 70, 300, 28, "Input Mesh Parameters :");
+          mesh_para_input->callback((Fl_Callback*)cb_mesh_para_input);
+          mesh_para_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Input* mesh_para_input
+        { mesh_file_btn = new Fl_Button(360, 70, 100, 28, "Mesh File");
+          mesh_file_btn->callback((Fl_Callback*)cb_mesh_file_btn);
+        } // Fl_Button* mesh_file_btn
+        { mod_para_file_btn = new Fl_Button(170, 160, 160, 28, "Add Model From File");
+          mod_para_file_btn->callback((Fl_Callback*)cb_mod_para_file_btn);
+        } // Fl_Button* mod_para_file_btn
+        { mod_ele_input_build = new Fl_Input(40, 445, 215, 28, "Input Model Element Data Name:");
+          mod_ele_input_build->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Input* mod_ele_input_build
+        { build_mod_btn = new Fl_Button(40, 545, 420, 28, "Build Model !");
+          build_mod_btn->callback((Fl_Callback*)cb_build_mod_btn);
+        } // Fl_Button* build_mod_btn
+        { mesh_para_output = new Fl_Output(40, 120, 300, 28, "Mesh Parameters :");
+          mesh_para_output->value("Unset");
+          mesh_para_output->box(FL_FLAT_BOX);
+          mesh_para_output->color(FL_BACKGROUND_COLOR);
+          mesh_para_output->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Output* mesh_para_output
+        { add_mod_btn = new Fl_Button(40, 160, 100, 28, "Add Model");
+          add_mod_btn->callback((Fl_Callback*)cb_add_mod_btn);
+        } // Fl_Button* add_mod_btn
+        { del_mod_btn = new Fl_Button(360, 160, 100, 28, "Delete Model");
+          del_mod_btn->callback((Fl_Callback*)cb_del_mod_btn);
+          del_mod_btn->deactivate();
+        } // Fl_Button* del_mod_btn
+        { rm_emp_bok_check = new Fl_Check_Button(280, 447, 180, 28, "Remove Empty Blocks");
+          rm_emp_bok_check->down_box(FL_DOWN_BOX);
+          rm_emp_bok_check->callback((Fl_Callback*)cb_rm_emp_bok_check);
+        } // Fl_Check_Button* rm_emp_bok_check
+        { mod_para_brw = new Fl_Browser(40, 205, 420, 210);
+          mod_para_brw->callback((Fl_Callback*)cb_mod_para_brw);
+          mod_para_brw->type(FL_MULTI_BROWSER);
+        } // Fl_Browser* mod_para_brw
+        { mod_file_out_btn = new Fl_Button(360, 500, 100, 28, "Model File");
+          mod_file_out_btn->callback((Fl_Callback*)cb_mod_file_out_btn);
+        } // Fl_Button* mod_file_out_btn
+        { mod_out_file_output = new Fl_Output(40, 500, 300, 28, "Output File Name :");
+          mod_out_file_output->box(FL_FLAT_BOX);
+          mod_out_file_output->color(FL_BACKGROUND_COLOR);
+          mod_out_file_output->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Output* mod_out_file_output
+        model_tab->end();
+      } // Fl_Group* model_tab
+      { forward_tab = new Fl_Group(10, 40, 480, 550, "Forward Modeling");
+        forward_tab->hide();
+        { grav_group = new Fl_Group(40, 340, 190, 58, "Forward Gravitational Data :");
+          grav_group->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+          { Vz_check = new Fl_Check_Button(40, 340, 50, 28, "Vz");
+            Vz_check->down_box(FL_DOWN_BOX);
+          } // Fl_Check_Button* Vz_check
+          { Vzx_check = new Fl_Check_Button(105, 340, 50, 28, "Vzx");
+            Vzx_check->down_box(FL_DOWN_BOX);
+          } // Fl_Check_Button* Vzx_check
+          { Vzy_check = new Fl_Check_Button(170, 340, 50, 28, "Vzy");
+            Vzy_check->down_box(FL_DOWN_BOX);
+          } // Fl_Check_Button* Vzy_check
+          { Vzz_check = new Fl_Check_Button(40, 370, 50, 28, "Vzz");
+            Vzz_check->down_box(FL_DOWN_BOX);
+          } // Fl_Check_Button* Vzz_check
+          grav_group->end();
+        } // Fl_Group* grav_group
+        { mod_file_btn = new Fl_Button(360, 70, 100, 28, "Model File");
+          mod_file_btn->callback((Fl_Callback*)cb_mod_file_btn);
+        } // Fl_Button* mod_file_btn
+        { obs_file_btn = new Fl_Button(360, 130, 100, 28, "Observe File");
+          obs_file_btn->callback((Fl_Callback*)cb_obs_file_btn);
+        } // Fl_Button* obs_file_btn
+        { mod_file_input = new Fl_Input(40, 70, 300, 28, "Input Model FIle :");
+          mod_file_input->callback((Fl_Callback*)cb_mod_file_input);
+          mod_file_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Input* mod_file_input
+        { mod_file_output = new Fl_Output(40, 240, 200, 28, "Chosen Model File :");
+          mod_file_output->box(FL_FLAT_BOX);
+          mod_file_output->color(FL_BACKGROUND_COLOR);
+          mod_file_output->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+          mod_file_output->value("Unset");
+        } // Fl_Output* mod_file_output
+        { obs_file_input = new Fl_Input(40, 130, 300, 28, "Input Observe  FIle :");
+          obs_file_input->callback((Fl_Callback*)cb_obs_file_input);
+          obs_file_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+          obs_file_input->tooltip("<xmin>/<dx>/<xmax>/<ymin>/<dy>/<ymax>/<elevation>|<filename>");
+        } // Fl_Input* obs_file_input
+        { obs_file_output = new Fl_Output(260, 240, 200, 28, "Chosen Observe File :");
+          obs_file_output->box(FL_FLAT_BOX);
+          obs_file_output->color(FL_BACKGROUND_COLOR);
+          obs_file_output->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+          obs_file_output->value("Unset");
+        } // Fl_Output* obs_file_output
+        { mod_ele_input = new Fl_Input(260, 315, 200, 28, "Model Element Data Name :");
+          mod_ele_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Input* mod_ele_input
+        { mag_group = new Fl_Group(40, 440, 420, 90);
+          { DeltaT_check = new Fl_Check_Button(40, 475, 65, 25, "DeltaT");
+            DeltaT_check->down_box(FL_DOWN_BOX);
+            DeltaT_check->deactivate();
+          } // Fl_Check_Button* DeltaT_check
+          { DeltaTx_check = new Fl_Check_Button(120, 475, 70, 25, "DeltaTx");
+            DeltaTx_check->down_box(FL_DOWN_BOX);
+            DeltaTx_check->deactivate();
+          } // Fl_Check_Button* DeltaTx_check
+          { DeltaTy_check = new Fl_Check_Button(200, 475, 70, 25, "DeltaTy");
+            DeltaTy_check->down_box(FL_DOWN_BOX);
+            DeltaTy_check->deactivate();
+          } // Fl_Check_Button* DeltaTy_check
+          { DeltaTz_check = new Fl_Check_Button(280, 475, 70, 25, "DeltaTz");
+            DeltaTz_check->down_box(FL_DOWN_BOX);
+            DeltaTz_check->deactivate();
+          } // Fl_Check_Button* DeltaTz_check
+          { Hax_check = new Fl_Check_Button(40, 505, 65, 25, "Hax");
+            Hax_check->down_box(FL_DOWN_BOX);
+            Hax_check->deactivate();
+          } // Fl_Check_Button* Hax_check
+          { Hay_check = new Fl_Check_Button(120, 505, 65, 25, "Hay");
+            Hay_check->down_box(FL_DOWN_BOX);
+            Hay_check->deactivate();
+          } // Fl_Check_Button* Hay_check
+          { Za_check = new Fl_Check_Button(200, 505, 65, 25, "Za");
+            Za_check->down_box(FL_DOWN_BOX);
+            Za_check->deactivate();
+          } // Fl_Check_Button* Za_check
+          { mag_para_input = new Fl_Input(220, 440, 240, 28, "Magnetization Parameters :  ");
+            mag_para_input->tooltip("<geo-inclination>/<geo-declination>/<mag-inclination>/<mag-declination>");
+            mag_para_input->deactivate();
+          } // Fl_Input* mag_para_input
+          mag_group->end();
+        } // Fl_Group* mag_group
+        { mag_data_check = new Fl_Check_Button(40, 405, 170, 30, "Forward Magnetic Data");
+          mag_data_check->down_box(FL_DOWN_BOX);
+          mag_data_check->callback((Fl_Callback*)cb_mag_data_check);
+        } // Fl_Check_Button* mag_data_check
+        { cal_btn = new Fl_Button(40, 545, 420, 28, "Calculate !");
+          cal_btn->callback((Fl_Callback*)cb_cal_btn);
+        } // Fl_Button* cal_btn
+        { noise_para_input = new Fl_Input(260, 370, 200, 28);
+          noise_para_input->tooltip("<noise-mean>/<noise-deviation>");
+          noise_para_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+          noise_para_input->deactivate();
+        } // Fl_Input* noise_para_input
+        { noise_check = new Fl_Check_Button(260, 345, 180, 30, "Input Noise Parameters :");
+          noise_check->down_box(FL_DOWN_BOX);
+          noise_check->callback((Fl_Callback*)cb_noise_check);
+        } // Fl_Check_Button* noise_check
+        { res_file_input = new Fl_Input(40, 190, 300, 28, "Input Prefix of Output FIle :");
+          res_file_input->callback((Fl_Callback*)cb_res_file_input);
+          res_file_input->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+        } // Fl_Input* res_file_input
+        { res_file_btn = new Fl_Button(360, 190, 100, 28, "Result File");
+          res_file_btn->callback((Fl_Callback*)cb_res_file_btn);
+        } // Fl_Button* res_file_btn
+        { res_file_output = new Fl_Output(40, 290, 200, 28, "Prefix of Output File :");
+          res_file_output->box(FL_FLAT_BOX);
+          res_file_output->color(FL_BACKGROUND_COLOR);
+          res_file_output->align(Fl_Align(FL_ALIGN_TOP_LEFT));
+          res_file_output->value("Unset");
+        } // Fl_Output* res_out_file_output
+        forward_tab->end();
+      } // Fl_Group* forward_tab
+      main_tabs->end();
+    } // Fl_Tabs* main_tabs
+    main_window->end();
+  } // Fl_Double_Window* main_window
+  main_window->show(argc, argv);
+  return Fl::run();
+}

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;
+}

test/right_mesh_para.txt

@@ -0,0 +1,3 @@
+10 20 990
+10 20 990
+10 20 490

test/right_mod_para.txt

@@ -0,0 +1 @@
+regular_block add 1.0 400/600/400/600/200/400

test/wrong_mesh_para.txt

@@ -0,0 +1,3 @@
+10 20 990
+10 20 990
+10 20

test/wrong_mod_para.txt

@@ -0,0 +1,2 @@
+this is wrong
+regular_block add 1.0 400/600/400/600/200/400