Merge pull request 'dev_yi' (#3) from dev_yi into main

Reviewed-on: #3

CMakeLists.txt

@@ -1,15 +1,9 @@
 cmake_minimum_required(VERSION 3.15.2)
 # 设置项目名称与语言
-project(GCTL VERSION 2.0)
+project(GCTL VERSION 1.0)
 # 添加配置配件编写的函数
 include(CMakePackageConfigHelpers)
 
-set(CMAKE_CXX_STANDARD 17)
-set(CMAKE_CXX_STANDARD_REQUIRED ON)
-
-# ExprTK库在macOS 15.4中编译会有错误 添加以下命令对应检查项
-add_compile_options(-Wno-missing-template-arg-list-after-template-kw)
-
 # 添加编译选项
 option(GCTL_OPENMP "Use the OpenMP library" ON)
 option(GCTL_NETCDF "Use the NetCDF library" ON)
@@ -17,7 +11,6 @@ option(GCTL_FFTW3 "Use the FFTW3 library" ON)
 option(GCTL_EEMD "Use the EEMD library" ON)
 option(GCTL_OPENBLAS "Use the Openblas library" OFF)
 option(GCTL_EXPRTK "Use the ExprTK library" ON)
-option(GCTL_GMT "Use the GMT library" ON)
 option(GCTL_CHECK_BOUNDER "Check array's index" OFF)
 option(GCTL_CHECK_SIZE "Check array's size" OFF)
 # 传递安装地址给编译期宏变量
@@ -33,26 +26,17 @@ message(STATUS "[GCTL] Use the FFTW3 library: " ${GCTL_FFTW3})
 message(STATUS "[GCTL] Use the EEMD library: " ${GCTL_EEMD})
 message(STATUS "[GCTL] Use the Openblas library: " ${GCTL_OPENBLAS})
 message(STATUS "[GCTL] Use the ExprTK library: " ${GCTL_EXPRTK})
-message(STATUS "[GCTL] Use the GMT library: " ${GCTL_GMT})
 message(STATUS "[GCTL] Check Bounder: " ${GCTL_CHECK_BOUNDER})
 message(STATUS "[GCTL] Check Size: " ${GCTL_CHECK_SIZE})
 
-find_library(NCURSES_LIB ncurses REQUIRED)
-
 if(GCTL_FFTW3)
 	if(NOT FFTW3_FOUND)
 		find_package(FFTW3 REQUIRED)
 		message(STATUS "Found FFTW3")
-		include_directories(${FFTW3_INCLUDE_DIRS})
+		include_directories(${FFTW3_INC_DIR})
 	endif()
 endif()
 
-if(GCTL_GMT)
-	find_package(GMT REQUIRED)
-	message(STATUS "Found GMT")
-	include_directories(${GMT_INC_DIR})
-endif()
-
 if(GCTL_EEMD)
 	if(NOT EEMD_FOUND)
 		find_package(EEMD REQUIRED)
@@ -70,8 +54,6 @@ if(GCTL_OPENBLAS)
 endif()
 
 if(GCTL_NETCDF)
-	find_package(HDF5)
-
 	if(NOT netCDF_FOUND)
 		find_package(netCDF REQUIRED)
 		include_directories(${netCDF_INCLUDE_DIR})

GCTLConfig.cmake.in

@@ -16,7 +16,6 @@ set(@PROJECT_NAME@_EEMD @GCTL_EEMD@)
 set(@PROJECT_NAME@_OPENMP @GCTL_OPENMP@)
 set(@PROJECT_NAME@_OPENBLAS @GCTL_OPENBLAS@)
 set(@PROJECT_NAME@_EXPRTK @GCTL_EXPRTK@)
-set(@PROJECT_NAME@_GMT @GCTL_GMT@)
 set(@PROJECT_NAME@_CHECK_BOUNDER @GCTL_CHECK_BOUNDER@)
 set(@PROJECT_NAME@_CHECK_SIZE @GCTL_CHECK_SIZE@)
 
@@ -26,15 +25,12 @@ message(STATUS "[GCTL] Use the EEMD library: " @GCTL_EEMD@)
 message(STATUS "[GCTL] Use the OpenMP library: " @GCTL_OPENMP@)
 message(STATUS "[GCTL] Use the Openblas library: " @GCTL_OPENBLAS@)
 message(STATUS "[GCTL] Use the ExprTK library: " @GCTL_EXPRTK@)
-message(STATUS "[GCTL] Use the GMT library: " @GCTL_GMT@)
 message(STATUS "[GCTL] Check Bounder: " @GCTL_CHECK_BOUNDER@)
 message(STATUS "[GCTL] Check Size: " @GCTL_CHECK_SIZE@)
 
 set(CMAKE_CXX_STANDARD 11)
 
 if(@PROJECT_NAME@_NETCDF)
-    find_package(HDF5)
-
     if(NOT netCDF_FOUND)
         find_package(netCDF REQUIRED)
         include_directories(${netCDF_INCLUDE_DIR})
@@ -49,17 +45,10 @@ endif()
 if(@PROJECT_NAME@_FFTW3)
     if(NOT FFTW3_FOUND)
         find_package(FFTW3 REQUIRED)
-        include_directories(${FFTW3_INCLUDE_DIRS})
+        include_directories(${FFTW3_INC_DIR})
     endif()
 endif()
 
-if(@PROJECT_NAME@_GMT)
-	if(NOT GMT_FOUND)
-		find_package(GMT REQUIRED)
-		include_directories(${GMT_INC_DIR})
-	endif()
-endif()
-
 if(@PROJECT_NAME@_EEMD)
     if(NOT EEMD_FOUND)
         find_package(EEMD REQUIRED)

README.md

@@ -2,8 +2,6 @@
 
 # Geophysical Computational Tools & Library (GCTL)
 
-> **注意**: 本文档由 Cursor AI 自动生成。文档中的函数名称、参数和用法说明可能存在错误，请以头文件（`.h`）中的实际函数声明为准。如发现任何不一致，请以源代码为准。
-
 GCTL 是一个用于地球物理研究的计算工具和 C++ 库。完整的软件包由核心库和额外的库以及命令行工具组成。本库采用现代 C++ 设计，提供高性能的数值计算和数据处理功能。
 
 ## 主要特性

example/CMakeLists.txt

@@ -23,17 +23,10 @@ add_example(windowfunc_ex OFF)
 add_example(legendre_ex OFF)
 add_example(refellipsoid_ex OFF)
 add_example(kde_ex OFF)
-add_example(meshio_ex OFF)
+add_example(meshio_ex ON)
 add_example(autodiff_ex OFF)
 add_example(multinary_ex OFF)
 add_example(text_io_ex OFF)
 add_example(getoption_ex OFF)
 add_example(process_ex OFF)
-add_example(array_ex OFF)
+add_example(array_ex OFF)
-add_example(gmt_ex OFF)
-add_example(gnuplot_ex OFF)
-add_example(cliplot_ex OFF)
-add_example(stl_io_ex OFF)
-add_example(ply_io_ex OFF)
-add_example(sparray_ex OFF)
-add_example(sparray2d_ex OFF)

example/cliplot_ex.cpp

@@ -1,47 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#include "../lib/graphic/cliplot.h"
-
-int main(int argc, char *argv[])
-{
-    double xmin = -3;
-    double xmax = 3;
-    double ymin = -1;
-    double ymax = 1;
-
-    gctl::cliplot c(81, 16, xmin, xmax, ymin, ymax);
-
-    //c.set_new_screen(true);
-    c.set_axis(5, 5);
-    c.set_digs(4);
-    c.set_wname("Time (s)");
-    c.set_hname("Value");
-    c.plot_func([](double x)->double{return sin(x);}, '.', GCTL_BOLDRED);
-    c.display();
-    return 0;
-}

example/eemd_ex.cpp

@@ -45,7 +45,7 @@
 
 #include "../lib/core.h"
 #include "../lib/io.h"
-#include "../lib/algorithms.h"
+#include "../lib/algorithm.h"
 
 using namespace gctl;
 

example/gmt_ex.cpp

@@ -1,65 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#include "../lib/graphic/gmt.h"
-
-int main(int argc, char *argv[])
-{
-    gctl::gmt gt;
-    /*
-    gt.begin_session("line_plot", "eps,png");
-    gt.call_module("basemap", "-R0/10/0/10 -JX5i/5i -Baf");
-    gt.call_module("plot", "tmp/data.txt -R0/10/0/10 -JX5i/5i -W1p,blue");
-    gt.end_session(false);
-    gt.save_session("line_plot");
-    */
-
-    // example 2
-    gctl::array<double> data(51*41);
-
-	double dist;
-	for (int i = 0; i < 51; ++i)
-	{
-		for (int j = 0; j < 41; ++j)
-		{
-			dist = sqrt((j-20)*(j-20) + (i-25)*(i-25));
-			data[j+i*41] = sin(dist/GCTL_Pi);
-		}
-	}
-
-    gt.begin_session("image_plot", "eps");
-    gt.call_module("set", "FONT_ANNOT_PRIMARY=10.5p,Times-Roman,black");
-
-    std::string gname = gt.create_grid(data, 41, 51, 0.0, 1.0, 0.0, 1.0);
-
-    gt.call_module("grd2cpt", gname + " -R0/40/0/50 -Crainbow -Z -D");
-    gt.call_module("grdimage", gname + " -R0/40/0/50 -JX1.5i/1.5i -X0.5i -Y0.5i -Bxag+l\"x (m)\" -Byag+l\"y (m)\"");
-    gt.call_module("psscale", "-Bxa -By+lm -Dx0.1i/-0.2i+w1.3i/0.05i+h");
-    gt.end_session();
-    gt.save_session("image_plot");
-    return 0;
-}

example/gnuplot_ex.cpp

@@ -1,97 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#include "../lib/core/macro.h"
-#include "../lib/graphic/gnuplot.h"
-
-int main(int argc, char *argv[])
-{
-    gctl::gnuplot gt;
-
-    //one line test
-    //gt.send("set terminal dumb");
-    //gt.send("plot [-pi/2:pi] cos(x),-(sin(x) > sin(x+1) ? sin(x) : sin(x+1))");
-
-    //buffer test (bessel animation)
-    /*
-    gt.to_buffer();
-    gt.send("set terminal gif animate delay 10 size 600,400");
-    gt.send("set output 'bessel.gif'");
-    gt.send("set palette rgb 3,9,9");
-    gt.send("unset key; unset colorbox; unset border; unset tics");
-    gt.send("set lmargin at screen 0.03");
-    gt.send("set bmargin at screen 0");
-    gt.send("set rmargin at screen 0.97");
-    gt.send("set tmargin at screen 1");
-    gt.send("set parametric", true);
-    gt.send("bessel(x,t) = besj0(x) * cos(2*pi*t)");
-    gt.send("n = 6 # number of zeros");
-    gt.send("k = (n*pi-1.0/4*pi)");
-    gt.send("u_0 = k + 1/(8*k) - 31/(384*k)**3 + 3779/(15360*k)**5");
-    gt.send("set urange [0:u_0]");
-    gt.send("set vrange[0:1.5*pi]");
-    gt.send("set cbrange [-1:1]");
-    gt.send("set zrange[-1:1]");
-    gt.send("set isosamples 200,100");
-    gt.send("set pm3d depthorder");
-    gt.send("set view 40,200");
-
-    std::string cmd;
-    for (float t = 0.0f; t < 2.0f; t += 0.02f)
-    {
-        cmd = "splot u*sin(v),u*cos(v),bessel(u," + std::to_string(t) + ") w pm3d ls 1";
-        gt.send(cmd);
-    }
-    gt.send("set output");
-
-    gt.save_buffer("bessel");
-    gt.send_buffer();
-    */
-
-    //data test
-    std::vector<double> x(101);
-    std::vector<double> y1(101);
-    std::vector<double> y2(101);
-    for (size_t i = 0; i < 101; i++)
-    {
-        x[i] = -1.0*GCTL_Pi + 2.0*GCTL_Pi*i/100.0;
-        y1[i] = sin(x[i]);
-        y2[i] = cos(x[i]);
-    }
-
-    std::vector<std::vector<double> > data;
-    data.push_back(x);
-    data.push_back(y1);
-    data.push_back(y2);
-
-    gt.to_buffer();
-    gt.send_data("$Data", data);
-    gt.send("plot $Data using 1:2 with lines title 'y1', '' using 1:3 with lines title 'y2'");
-    gt.save_buffer("sin_cos");
-    gt.send_buffer();
-    return 0;
-}

example/meshio_ex.cpp

@@ -0,0 +1,87 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "../lib/core.h"
+#include "../lib/io.h"
+
+using namespace gctl;
+
+int main(int argc, char const *argv[]) try
+{
+    mesh_io mshio;
+
+/*  
+    mshio.read_tetgen_ascii("tmp/ex1.1");
+    mshio.edit_group(Disable, GeometryTag, 5);
+    mshio.edit_group(GeometryTag, 1, PhysicalTag, 1);
+    mshio.edit_group(GeometryTag, 2, PhysicalTag, 2);
+    mshio.edit_group(GeometryTag, 3, PhysicalTag, 3);
+    mshio.edit_group(GeometryTag, 4, PhysicalTag, 4);
+    mshio.edit_group(GeometryTag, 1, "Boundary");
+    mshio.edit_group(GeometryTag, 2, "Body1");
+    mshio.edit_group(GeometryTag, 3, "Body2");
+    mshio.edit_group(GeometryTag, 4, "Body3");
+    mshio.save_gmsh_v2_ascii("tmp/ex1.1");
+*/
+
+    mshio.read_gmsh_v2_ascii("tmp/ex1.1");
+    mshio.convert_tags_to_data(GeometryTag);
+
+    array<double> body_val(mshio.element_size("Body2"), 2.0);
+    mshio.add_element_data("BodyValue", "Body2", body_val);
+
+    array<double> body_val2(mshio.element_size("Body3"), 1.0);
+    mshio.add_element_data("BodyValue", "Body3", body_val2);
+
+    mshio.save_gmsh_v2_ascii("tmp/ex1.2");
+    //mshio.save_vtk_legacy_ascii("tmp/ex1.1");
+    mshio.info();
+
+    const array<vertex3dc> &nodes = mshio.get_nodes();
+
+    array<tetrahedron> body2_tets;
+    mshio.export_elements_to(body2_tets, "All");
+
+    gmshio gio;
+    gio.init_file("tmp.msh", Output);
+    gio.set_packed(NotPacked, Output);
+    gio.save_mesh(body2_tets, nodes);
+/*
+    mshio.read_gmsh_v2_ascii("tmp/wjb.1");
+    mshio.edit_group(Disable);
+    mshio.edit_group(Enable, GeometryTag, 3);
+    mshio.edit_group(Enable, GeometryTag, 8);
+    mshio.edit_group(Enable, GeometryTag, 9);
+
+    mshio.save_gmsh_v2_ascii("tmp/wjb.2");
+*/
+    return 0;
+}
+catch(std::exception &e)
+{
+    GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}

example/ply_io_ex.cpp

@@ -1,44 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
- #include "../lib/core.h"
- #include "../lib/io.h"
- 
- using namespace gctl;
- 
- int main(int argc, char const *argv[]) try
- {
-     array<vertex3dc> nodes;
-     array<triangle> tris;
-     read_ply_binary("tmp_binary", nodes, tris);
-     save_ply_ascii("tmp", nodes, tris);
-     return 0;
- }
- catch(std::exception &e)
- {
-     GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
- }

example/refellipsoid_ex.cpp

@@ -25,48 +25,20 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#include "../lib/geometry.h"
+#include "../lib/geodesy.h"
 
 using namespace gctl;
 
 int main(int argc, char const *argv[]) try
 {
-    refellipsoid ellip;
+    refellipsoid relli(WGS84);
-    ellip.set(WGS84);
+    refellipsoid relli2(1000, 900, false);
-/*
+    refellipsoid relli3(1000, 10, true);
-    point3ds ps;
-    ps.set_io_precision(17);
-    ps.lon = -110;
-    ellip.geodetic2spherical(-40, 0, ps.lat, ps.rad);
-    std::cout << "Geocentric: " << ps << std::endl;
 
-    point3dc pc = ps.s2c();
+    for (size_t i = 0; i <= 90; i++)
-    pc.set_io_precision(17);
+    {
-    std::cout << "XYZ: " <<  pc << std::endl;
+        std::cout << i << " " << relli.radis_at(1.0*i) << " " << relli2.radis_at(1.0*i) << " " << relli3.radis_at(1.0*i) << "\n";
-
+    }
-    ellip.spherical2geodetic(ps, ps.lon, ps.lat, ps.rad);
-    std::cout << "Geodetic: " << ps << std::endl;
-
-    ellip.xyz2geodetic(pc, ps.lon, ps.lat, ps.rad);
-    std::cout << "Geodetic: " << ps << std::endl;
-*/
-
-    point3ds ps;
-    ps.set_io_precision(17);
-    // 大地坐标 (110,40,0)
-    ps.lon = 110;
-    ellip.geodetic2spherical(40.0, 0.0, ps.lat, ps.rad);
-    // 球心坐标 (110,39.810615323061491,6369344.5441424493)
-    std::cout << "Geocentric: " << ps << std::endl;
-    // 转换为大地坐标
-    point3ds pd;
-    ellip.xyz2geodetic(ps.s2c(), pd.lon, pd.lat, pd.rad);
-    std::cout << "Geodetic: " << pd << std::endl;
-    // 400km高
-    pd.rad = 400000.0;
-    ellip.geodetic2spherical(pd.lat, pd.rad, pd.lat, pd.rad);
-    std::cout << "Geocentric: " << pd << std::endl;
-    return 0;
 }
 catch(std::exception &e)
 {

example/sparray_ex.cpp

@@ -1,114 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#include "../lib/core.h"
-#include "../lib/algorithms.h"
-
-int main(int argc, char const *argv[])
-{
-	srand(time(0));
-	gctl::sparray<double> M(50, 0.0);
-	gctl::sparray<double> L;
-
-	int tmp_id[13] = {4, 2, 6, 7, 8, 12, 23, 43, 33, 47, 38, 15, 1};
-	int tmp_size = 13;
-	for (int i = 0; i < tmp_size; i++)
-	{
-		M.set(tmp_id[i], gctl::random(10.0, 20.0));
-	}
-
-	M.show_list();
-	std::cout << "************" << std::endl;
-
-	M.copy_to(L, 2);
-	L.show_list();
-
-	M.set(18, 100);
-	M.set(38, 100);
-	M.remove(12);
-
-	std::cout << "************" << std::endl;
-	M.show_list();
-
-	gctl::array<double> C(50, -1.0);
-	M.export_dense(C, 1.0, gctl::AppendVal);
-	for (int i = 0; i < C.size(); i++)
-	{
-		std::cout << C.at(i) << " ";
-	}
-	std::cout << std::endl;
-
-	gctl::sparray<double> N(C, -1, 1e-10);
-	N.show_list();
-
-	std::cout << "Test 2" << std::endl;
-
-	M.clear();
-	M.malloc(20, 0.0);
-
-	M.set(5, 6.5);
-	M.set(9, 9.1);
-	M.set(3, 4.3);
-	M.set(17, 1.4);
-	M.set(10, 3.5);
-	M.set(7, 7.4);
-
-	M.show_list();
-
-	M.remove(17);
-	M.remove(9);
-
-	std::cout << "**********" << std::endl;
-	M.show_list();
-
-	std::cout << "Test 3" << std::endl;
-
-	int test_size = 500000;
-
-	M.clear();
-	M.malloc(test_size, 0.0);
-
-	clock_t start = clock();
-	for (int i = 0; i < 300000; i++)
-	{
-		M.set(i, gctl::random(10.0, 20.0));
-	}
-	for (int i = 300001; i < test_size; i++)
-	{
-		M.set(i, gctl::random(10.0, 20.0));
-	}
-	M.set(300000, 15.8888);
-
-	clock_t end = clock();
-	std::cout << "sparray set() time: " << 1000.0*(end - start)/(double)CLOCKS_PER_SEC << " ms" << std::endl;
-
-	std::cout << "M.at(300000) = " << M.value(300000) << std::endl;
-
-	M.show_list(299995, 300005);
-
-	return 0;
-}

example/text_io_ex.cpp

@@ -25,6 +25,7 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
+#include "../lib/core.h"
 #include "../lib/io.h"
 
 using namespace gctl;
@@ -32,20 +33,26 @@ using namespace gctl;
 int main(int argc, char const *argv[]) try
 {
 
-    dsv_io tc;
+    dsv_io tc, tout;
-    tc.delimeter('|');
+    tc.set_delimeter('|');
-    tc.head_number(1);
+    tc.load_text("tmp/world_data", ".txt", BothHead);
-    tc.load_text("tmp/world_data", ".txt", ColHead|RowHead);
+    //tc.info();
-    tc.info(AttInfo|HeadInfo|TagInfo);
 
-    tc.filter("America", "Continent_s", ColHead);
+    //tc.set_column_type(Int, "IndepYear_n");
-    //tc.filter("America", "BLZ", RowHead);
+    //tc.filt_column("IndepYear_n < 0", {"IndepYear_n"}, {"Name_s", "Population_n", "GNP_n"}, tout);
-    //tc.save_text("out");
 
-    dsv_io tc2 = tc.export_table();
+    tc.filt_column("America", "Continent_s", {"Name_s", "Population_n", "GNP_n"}, tout);
-    //tc2.head_records(tc.head_records());
+    //tc.match_column("America", "Continent_s", {}, tout);
-    tc2.delimeter('|');
+
-    tc2.save_text("out");
+    //tout.add_column("GNP_n2", "Population_n");
+    //array<int> GNP_n2(tout.row_number(), 1000.0);
+    //tout.fill_column(GNP_n2, "GNP_n2");
+
+    int lr_id = tout.add_row();
+    tout.fill_row(array<std::string>{"Asia", "China", "14000000", "1949"}, lr_id);
+
+    tout.set_delimeter('|');
+    tout.save_text("out");
 
 /*
     geodsv_io tc;

lib/CMakeLists.txt

@@ -5,11 +5,12 @@ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3")
 set(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib)
 
 # 设定库源文件文件夹
-aux_source_directory(core/ GCTL_CORE_SRC)
+aux_source_directory(geometry/ GCTL_GEOMETRY_SRC)
+aux_source_directory(algorithm/ GCTL_ALGORITHM_SRC)
 aux_source_directory(io/ GCTL_IO_SRC)
 aux_source_directory(utility/ GCTL_UTILITY_SRC)
-aux_source_directory(math/ GCTL_MATH_SRC)
+aux_source_directory(maths/ GCTL_MATHS_SRC)
-aux_source_directory(graphic/ GCTL_GRAPHIC_SRC)
+aux_source_directory(geodesy/ GCTL_GEODESY_SRC)
 # 在windows系统下使用改编过的getopt_win源文件替代GNU getopt
 # 若在linux系统下则不适用此源文件
 if(NOT WIN32)
@@ -23,11 +24,11 @@ endif()
 # 以下部分为库的编译
 # 注意目标名必须唯一 所以不能直接生成相同名称的动态库与静态库
 # 注意此处不必为目标名称添加lib前缀和相应后缀，cmake会自行添加
-add_library(gctl SHARED ${GCTL_CORE_SRC} ${GCTL_UTILITY_SRC} 
+add_library(gctl SHARED ${GCTL_UTILITY_SRC} ${GCTL_ALGORITHM_SRC} ${GCTL_GEODESY_SRC}
-	${GCTL_GRAPHIC_SRC} ${GCTL_IO_SRC} ${GCTL_MATH_SRC})
+	${GCTL_GEOMETRY_SRC} ${GCTL_IO_SRC} ${GCTL_MATHS_SRC})
 # 首先添加静态库的生成命令
-add_library(gctl_static STATIC ${GCTL_CORE_SRC} ${GCTL_UTILITY_SRC} 
+add_library(gctl_static STATIC ${GCTL_UTILITY_SRC} ${GCTL_ALGORITHM_SRC} ${GCTL_GEODESY_SRC}
-	${GCTL_GRAPHIC_SRC} ${GCTL_IO_SRC} ${GCTL_MATH_SRC})
+	${GCTL_GEOMETRY_SRC} ${GCTL_IO_SRC} ${GCTL_MATHS_SRC})
 # 设置静态库的输出名称从而获得与动态库名称相同的静态库
 set_target_properties(gctl_static PROPERTIES OUTPUT_NAME "gctl")
 # 设置输出目标属性以同时输出动态库与静态库
@@ -41,29 +42,19 @@ set_target_properties(gctl_static PROPERTIES INSTALL_RPATH /usr/local/lib)
 set_target_properties(gctl PROPERTIES CXX_STANDARD 17 CXX_STANDARD_REQUIRED ON)
 set_target_properties(gctl_static PROPERTIES CXX_STANDARD 17 CXX_STANDARD_REQUIRED ON)
 
-target_link_libraries(gctl PUBLIC ${NCURSES_LIB})
-target_link_libraries(gctl_static ${NCURSES_LIB})
-
 if(GCTL_NETCDF)
 	target_link_libraries(gctl PUBLIC ${netCDF_CXX_LEGACY_LIB})
 	target_link_libraries(gctl_static ${netCDF_CXX_LEGACY_LIB})
 	target_link_libraries(gctl PUBLIC ${netCDF_LIBRARIES})
 	target_link_libraries(gctl_static ${netCDF_LIBRARIES})
-	target_link_libraries(gctl PUBLIC ${HDF5_LIBRARIES})
-	target_link_libraries(gctl_static ${HDF5_LIBRARIES})
 endif()
 
 # 连接动态库与静态库
 if(GCTL_FFTW3)
-	find_library(FFTW_LIB ${FFTW3_LIBRARIES} HINTS ${FFTW3_LIBRARY_DIRS})
+	# 因为是手动设置的动态库地址 使用find_library更保险
-	target_link_libraries(gctl PUBLIC ${FFTW_LIB})
+	find_library(LIB_TAR ${FFTW3_LIB} HINTS ${FFTW3_LIB_DIR})
-	target_link_libraries(gctl_static ${FFTW_LIB})
+	target_link_libraries(gctl PUBLIC ${LIB_TAR})
-endif()
+	target_link_libraries(gctl_static ${LIB_TAR})
-
-if(GCTL_GMT)
-	find_library(GMT_LIB_TAR ${GMT_LIB} HINTS ${GMT_LIB_DIR})
-	target_link_libraries(gctl PUBLIC ${GMT_LIB_TAR})
-	target_link_libraries(gctl_static ${GMT_LIB_TAR})
 endif()
 
 if(GCTL_EEMD)
@@ -111,17 +102,19 @@ endif()
 file(GLOB GCTL_HEAD *.h)
 file(GLOB GCTL_CORE_HEAD core/*.h)
 file(GLOB GCTL_IO_HEAD io/*.h)
-file(GLOB GCTL_MATH_HEAD math/*.h)
+file(GLOB GCTL_ALGORITHM_HEAD algorithm/*.h)
-file(GLOB GCTL_POLY_HEAD poly/*.h)
+file(GLOB GCTL_MATHS_HEAD maths/*.h)
+file(GLOB GCTL_GEOMETRY_HEAD geometry/*.h)
 file(GLOB GCTL_UTILITY_HEAD utility/*.h)
-file(GLOB GCTL_GRAPHIC_HEAD graphic/*.h)
+file(GLOB GCTL_MESH_HEAD mesh/*.h)
+file(GLOB GCTL_GEODESY_HEAD geodesy/*.h)
 
 if(NOT WIN32)
 	list(REMOVE_ITEM GCTL_UTILITY_HEAD "${PROJECT_SOURCE_DIR}/lib/utility/getopt_win.h")
 endif()
 
 if(NOT GCTL_FFTW3)
-	list(REMOVE_ITEM GCTL_MATH_HEAD "${PROJECT_SOURCE_DIR}/lib/math/fft.h")
+	list(REMOVE_ITEM GCTL_ALGORITHM_HEAD "${PROJECT_SOURCE_DIR}/lib/algorithm/fft.h")
 endif()
 
 if(NOT GCTL_NETCDF)
@@ -131,7 +124,9 @@ endif()
 install(FILES ${GCTL_HEAD} DESTINATION include/gctl)
 install(FILES ${GCTL_CORE_HEAD} DESTINATION include/gctl/core)
 install(FILES ${GCTL_IO_HEAD} DESTINATION include/gctl/io)
-install(FILES ${GCTL_MATH_HEAD} DESTINATION include/gctl/math)
+install(FILES ${GCTL_MATHS_HEAD} DESTINATION include/gctl/maths)
-install(FILES ${GCTL_POLY_HEAD} DESTINATION include/gctl/poly)
+install(FILES ${GCTL_ALGORITHM_HEAD} DESTINATION include/gctl/algorithm)
+install(FILES ${GCTL_GEOMETRY_HEAD} DESTINATION include/gctl/geometry)
 install(FILES ${GCTL_UTILITY_HEAD} DESTINATION include/gctl/utility)
-install(FILES ${GCTL_GRAPHIC_HEAD} DESTINATION include/gctl/graphic)
+install(FILES ${GCTL_MESH_HEAD} DESTINATION include/gctl/mesh)
+install(FILES ${GCTL_GEODESY_HEAD} DESTINATION include/gctl/geodesy)

lib/algorithm.h

@@ -1,4 +1,4 @@
-/********************************************************
+/********************************************************
  *  ██████╗  ██████╗████████╗██╗
  * ██╔════╝ ██╔════╝╚══██╔══╝██║
  * ██║  ███╗██║        ██║   ██║
@@ -25,21 +25,24 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#include "../lib/core.h"
+#ifndef _GCTL_ALGORITHM_H
-#include "../lib/io.h"
+#define _GCTL_ALGORITHM_H
 
-using namespace gctl;
+#include "algorithm/algorithm_func.h"
+#include "algorithm/glni.h"
+#include "algorithm/interpolate.h"
+#include "algorithm/extrapolate.h"
+#include "algorithm/heap_sort.h"
+#include "algorithm/boxsort2d.h"
+#include "algorithm/boxsort_sph.h"
+#include "algorithm/variogram.h"
+#include "algorithm/fir_filter.h"
+#include "algorithm/space_filter.h"
+#include "algorithm/sinkhorn.h"
+#include "algorithm/kde.h"
+#include "algorithm/autodiff.h"
+#include "algorithm/multinary.h"
+#include "algorithm/emd.h"
+#include "algorithm/fft_filter.h"
 
-int main(int argc, char const *argv[]) try
+#endif // _GCTL_ALGORITHM_H
-{
-    array<vertex3dc> nodes;
-    array<triangle> tris;
-    read_stl_binary("tmp/Stanford_Bunny_Binary", nodes, tris);
-    save_stl_ascii("tmp", tris, "Stanford_Bunny");
-    save_stl_binary("tmp_binary", tris, "Stanford_Bunny");
-    return 0;
-}
-catch(std::exception &e)
-{
-    GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
-}

lib/algorithm/algorithm_func.cpp

@@ -1,4 +1,4 @@
-/********************************************************
+/********************************************************
  *  ██████╗  ██████╗████████╗██╗
  * ██╔════╝ ██╔════╝╚══██╔══╝██║
  * ██║  ███╗██║        ██║   ██║
@@ -25,483 +25,7 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#include "gmath.h"
+#include "algorithm_func.h"
-
-int gctl::random(int low, int hig)
-{
-	return rand()%(hig - low + 1) + low;
-}
-
-double gctl::random(double low, double hig)
-{
-	double f = (double) rand()/RAND_MAX;
-	return f*(hig-low) + low;
-}
-
-std::complex<double> gctl::random(std::complex<double> low, std::complex<double> hig)
-{
-	std::complex<double> c;
-	double r = (double) rand()/RAND_MAX;
-	double i = (double) rand()/RAND_MAX;
-	double rh = std::max(hig.real(), low.real());
-	double rl = std::min(hig.real(), low.real());
-	double ih = std::max(hig.imag(), low.imag());
-	double il = std::min(hig.imag(), low.imag());
-
-	c.real(r*(rh - rl) + rl);
-	c.imag(i*(ih - il) + il);
-	return c;
-}
-
-bool gctl::isequal(double f, double v, double eps)
-{
-	return fabs(f - v) < eps;
-}
-
-double gctl::sign(double a)
-{
-	if (a > GCTL_ZERO) return 1.0;
-	if (a < -1.0*GCTL_ZERO) return -1.0; 
-	return 0.0;
-}
-
-//利用二分法求一个正数的n次方根 注意输入值小于1的情况
-double gctl::sqrtn(double d,int n,double eps)
-{
-	double xmin,xmax,halfx;
-	if (d == 1)
-	{
-		return d;
-	}
-	else if (d > 1)
-	{
-		xmin = 1;
-		xmax = d;
-		halfx = 0.5*(xmin+xmax);
-		while (fabs(d - pow(halfx,n)) > eps)
-		{
-			if (pow(halfx,n) > d)
-			{
-				xmax = halfx;
-				halfx = 0.5*(xmin+xmax);
-			}
-			else
-			{
-				xmin = halfx;
-				halfx = 0.5*(xmin+xmax);
-			}
-		}
-	}
-	else
-	{
-		xmin = 0;
-		xmax = 1;
-		halfx = 0.5*(xmin+xmax);
-		while (fabs(d - pow(halfx,n)) > eps)
-		{
-			if (pow(halfx,n) > d)
-			{
-				xmax = halfx;
-				halfx = 0.5*(xmin+xmax);
-			}
-			else
-			{
-				xmin = halfx;
-				halfx = 0.5*(xmin+xmax);
-			}
-		}
-	}
-	return halfx;
-}
-
-double gctl::geographic_area(double lon1, double lon2, double lat1, double lat2, double R)
-{
-	return fabs(R*R*(arc(lon2) - arc(lon1))*(sind(lat2) - sind(lat1)));
-}
-
-double gctl::geographic_distance(double lon1, double lon2, double lat1, double lat2, double R)
-{
-	double n1 = arc(lon1), n2 = arc(lon2), t1 = arc(lat1), t2 = arc(lat2);
-	double a = sin(0.5*(t2 - t1)), b = sin(0.5*(n2 - n1));
-	return 2*R*asin(sqrt(a*a + cos(t1)*cos(t2)*b*b));
-}
-
-double gctl::ellipse_radius_2d(double x_len, double y_len, double arc, double x_arc)
-{
-	if (fabs(x_len - y_len) < 1e-16) // 就是个圆 直接加
-	{
-		return 0.5*(x_len + y_len);
-	}
-	
-	return sqrt(power2(x_len*cos(arc - x_arc)) + power2(y_len*sin(arc - x_arc)));
-}
-
-void gctl::ellipse_plus_elevation_2d(double x_len, double y_len, double arc, double elev,
-	double &out_arc, double &out_rad, double x_arc)
-{
-	if (fabs(x_len - y_len) < 1e-8) // 就是个圆 直接加
-	{
-		out_arc = arc;
-		out_rad = 0.5*(x_len + y_len) + elev;
-		return;
-	}
-	
-	if (fabs(arc) < 1e-8) // 弧度太小 直接加和
-	{
-		out_arc = arc;
-		out_rad = x_len + elev;
-		return;
-	}
-
-	if (fabs(fabs(arc) - 0.5*GCTL_Pi) < 1e-8) // 到了弧顶 直接加和
-	{
-		out_arc = arc;
-		out_rad = y_len + elev;
-		return;
-	}
-	
-	double ellip_rad = ellipse_radius_2d(x_len, y_len, arc, x_arc);
-	double alpha = atan(x_len*x_len*sin(arc)/(y_len*y_len*cos(arc)));
-	double sin_alpha = sin(alpha);
-	double lx = ellip_rad * sin(alpha - arc)/sin_alpha;
-	double lr = ellip_rad * sin(arc)/sin_alpha + elev;
-	out_rad = sqrt(lx*lx + lr*lr - 2.0*lx*lr*cos(GCTL_Pi - alpha));
-	out_arc = acos((lx*lx + out_rad*out_rad - lr*lr)/(2.0*out_rad*lx));
-	if (arc < 0.0) out_arc *= -1.0;
-
-	return;
-}
-
-double gctl::ellipsoid_radius(double x_len, double y_len, double z_len, double phi, double theta)
-{
-	return x_len*y_len*z_len/sqrt(pow(y_len*z_len*sin(theta)*cos(phi),2) + 
-		pow(x_len*z_len*sin(theta)*sin(phi),2) + pow(x_len*y_len*cos(theta),2));
-}
-
-// 使用牛顿迭代法计算一个矩阵的近似逆
-double gctl::newton_inverse(const _2d_matrix &in_mat, _2d_matrix &inverse_mat, 
-	double epsilon, int iter_times, bool initiate)
-{
-	if (in_mat.empty())
-		throw runtime_error("The input matrix is empty. Thrown by gctl::newton_inverse(...)");
-
-	if (in_mat.row_size() != in_mat.col_size())
-		throw logic_error("The input matrix is square. Thrown by gctl::newton_inverse(...)");
-
-	if (iter_times < 0)
-		throw invalid_argument("Invalid iteration times. Thrown by gctl::newton_inverse(...)");
-
-	int m_size = in_mat.row_size();
-
-	if (initiate)
-	{
-		// 初始化逆矩阵 使用输入矩阵的对角元素构建初始矩阵
-		if (!inverse_mat.empty()) inverse_mat.clear();
-		inverse_mat.resize(m_size, m_size, 0.0);
-		for (int i = 0; i < m_size; i++)
-		{
-			inverse_mat[i][i] = 1.0/in_mat[i][i];
-		}
-	}
-
-	if (inverse_mat.row_size() != m_size || inverse_mat.col_size() != m_size)
-		throw logic_error("Invalid output matrix size. From gctl::newton_inverse(...)");
-
-	// 迭代的收敛条件 || E - in_mat*inverse_mat ||_inf < 1
-	// 计算矩阵的无穷大范数
-	double maxi_mod = 0.0, mod, ele;
-	for (int i = 0; i < m_size; i++)
-	{
-		mod = 0.0;
-		for (int j = 0; j < m_size; j++)
-		{
-			ele = 0.0;
-			for (int k = 0; k < m_size; k++)
-			{
-				ele += in_mat[i][k] * inverse_mat[k][j];
-			}
-
-			if (i == j)
-			{
-				mod += GCTL_FABS(1.0 - ele);
-			}
-			else
-			{
-				mod += GCTL_FABS(ele);
-			}
-		}
-
-		maxi_mod = GCTL_MAX(maxi_mod, mod);
-	}
-
-	if (maxi_mod >= 1.0)
-	{
-		GCTL_ShowWhatError("The iteration may not converge. From gctl::newton_inverse(...)", 
-			GCTL_WARNING_ERROR, 0, 0, 0);
-	}
-
-	array<double> col_ax(m_size, 0.0); // A*X 的列向量
-	_2d_matrix tmp_mat(m_size, m_size, 0.0);
-	for (int t = 0; t < iter_times; t++)
-	{
-		if (maxi_mod <= epsilon)
-			break;
-
-		for (int i = 0; i < m_size; i++)
-		{
-			for (int j = 0; j < m_size; j++)
-			{
-				tmp_mat[i][j] = inverse_mat[i][j];
-			}
-		}
-
-		for (int n = 0; n < m_size; n++)
-		{
-			for (int i = 0; i < m_size; i++)
-			{
-				for (int j = 0; j < m_size; j++)
-				{
-					col_ax[j] = 0.0;
-					for (int k = 0; k < m_size; k++)
-					{
-						col_ax[j] += in_mat[j][k] * tmp_mat[k][i];
-					}
-					col_ax[j] *= -1.0;
-				}
-				col_ax[i] += 2.0;
-
-				ele = 0.0;
-				for (int j = 0; j < m_size; j++)
-				{
-					ele += tmp_mat[n][j] * col_ax[j];
-				}
-
-				inverse_mat[n][i] = ele;
-			}
-		}
-
-		maxi_mod = 0.0;
-		for (int i = 0; i < m_size; i++)
-		{
-			mod = 0.0;
-			for (int j = 0; j < m_size; j++)
-			{
-				ele = 0.0;
-				for (int k = 0; k < m_size; k++)
-				{
-					ele += in_mat[i][k] * inverse_mat[k][j];
-				}
-
-				if (i == j)
-				{
-					mod += GCTL_FABS(1.0 - ele);
-				}
-				else
-				{
-					mod += GCTL_FABS(ele);
-				}
-			}
-
-			maxi_mod = GCTL_MAX(maxi_mod, mod);
-		}
-	}
-
-	return maxi_mod;
-}
-
-// 使用牛顿迭代法计算一个矩阵的近似逆
-double gctl::newton_inverse(const spmat<double> &in_mat, _2d_matrix &inverse_mat, 
-	double epsilon, int iter_times, bool initiate)
-{
-	if (in_mat.empty())
-		throw runtime_error("The input matrix is empty. Thrown by gctl::newton_inverse(...)");
-
-	if (in_mat.row_size() != in_mat.col_size())
-		throw logic_error("The input matrix is square. Thrown by gctl::newton_inverse(...)");
-
-	if (iter_times < 0)
-		throw invalid_argument("Invalid iteration times. Thrown by gctl::newton_inverse(...)");
-
-	int m_size = in_mat.row_size();
-
-	if (initiate)
-	{
-		// 初始化逆矩阵 使用输入矩阵的对角元素构建初始矩阵
-		if (!inverse_mat.empty()) inverse_mat.clear();
-		inverse_mat.resize(m_size, m_size, 0.0);
-		for (int i = 0; i < m_size; i++)
-		{
-			inverse_mat[i][i] = 1.0/in_mat.at(i, i);
-		}
-	}
-
-	if (inverse_mat.row_size() != m_size || inverse_mat.col_size() != m_size)
-		throw logic_error("Invalid output matrix size. From gctl::newton_inverse(...)");
-
-	// 迭代的收敛条件 || E - in_mat*inverse_mat ||_inf < 1
-	// 计算矩阵的无穷大范数
-	double maxi_mod = 0.0, mod, ele;
-	array<double> tmp_arr(m_size, 0.0); // A*X 的列向量
-	for (int i = 0; i < m_size; i++)
-	{
-		mod = 0.0;
-		for (int j = 0; j < m_size; j++)
-		{
-			for (int k = 0; k < m_size; k++)
-			{
-				tmp_arr[k] = inverse_mat[k][j];
-			}
-
-			ele = in_mat.multiply_vector(tmp_arr, i);
-
-			if (i == j)
-			{
-				mod += GCTL_FABS(1.0 - ele);
-			}
-			else
-			{
-				mod += GCTL_FABS(ele);
-			}
-		}
-
-		maxi_mod = GCTL_MAX(maxi_mod, mod);
-	}
-
-	if (maxi_mod >= 1.0)
-	{
-		GCTL_ShowWhatError("The iteration may not converge. From gctl::newton_inverse(...)", 
-			GCTL_WARNING_ERROR, 0, 0, 0);
-	}
-
-	array<double> col_ax(m_size, 0.0);
-	_2d_matrix tmp_mat(m_size, m_size, 0.0);
-	for (int t = 0; t < iter_times; t++)
-	{
-		if (maxi_mod <= epsilon)
-			break;
-		//else std::cout << "epsilon = " << maxi_mod << std::endl;
-
-		for (int i = 0; i < m_size; i++)
-		{
-			for (int j = 0; j < m_size; j++)
-			{
-				tmp_mat[i][j] = inverse_mat[i][j];
-			}
-		}
-
-		for (int n = 0; n < m_size; n++)
-		{
-			for (int i = 0; i < m_size; i++)
-			{
-				for (int j = 0; j < m_size; j++)
-				{
-					for (int k = 0; k < m_size; k++)
-					{
-						tmp_arr[k] = tmp_mat[k][i];
-					}
-					col_ax[j] = in_mat.multiply_vector(tmp_arr, j);
-					col_ax[j] *= -1.0;
-				}
-				col_ax[i] += 2.0;
-
-				ele = 0.0;
-				for (int j = 0; j < m_size; j++)
-				{
-					ele += tmp_mat[n][j] * col_ax[j];
-				}
-
-				inverse_mat[n][i] = ele;
-			}
-		}
-
-		maxi_mod = 0.0;
-		for (int i = 0; i < m_size; i++)
-		{
-			mod = 0.0;
-			for (int j = 0; j < m_size; j++)
-			{
-				for (int k = 0; k < m_size; k++)
-				{
-					tmp_arr[k] = inverse_mat[k][j];
-				}
-
-				ele = in_mat.multiply_vector(tmp_arr, i);
-
-				if (i == j)
-				{
-					mod += GCTL_FABS(1.0 - ele);
-				}
-				else
-				{
-					mod += GCTL_FABS(ele);
-				}
-			}
-
-			maxi_mod = GCTL_MAX(maxi_mod, mod);
-		}
-	}
-
-	return maxi_mod;
-}
-
-void gctl::schmidt_orthogonal(const array<double> &a, array<double> &e, int a_s)
-{
-	if (a.empty())
-		throw runtime_error("The input array is empty. Thrown by gctl::schmidt_orthogonal(...)");
-
-	if (a_s <= 1) // a_s >= 2
-		throw invalid_argument("vector size must be bigger than one. Thrown by gctl::schmidt_orthogonal(...)");
-
-	int t_s = a.size();
-	if (t_s%a_s != 0 || t_s <= 3) // t_s >= 4
-		throw invalid_argument("incompatible total array size. Thrown by gctl::schmidt_orthogonal(...)");
-
-	int len = t_s/a_s;
-	if (len < a_s)
-		throw invalid_argument("the vectors are over-determined. Thrown by gctl::schmidt_orthogonal(...)");
-
-	e.resize(t_s, 0.0);
-
-	double ae, ee;
-	for (int i = 0; i < a_s; i++)
-	{
-		for (int l = 0; l < len; l++)
-		{
-			e[l + i*len] = a[l + i*len];
-		}
-
-		for (int m = 0; m < i; m++)
-		{
-			ae = ee = 0.0;
-			for (int n = 0; n < len; n++)
-			{
-				ae += a[n + i*len] * e[n + m*len];
-				ee += e[n + m*len] * e[n + m*len];
-			}
-
-			for (int n = 0; n < len; n++)
-			{
-				e[n + i*len] -= e[n + m*len] * ae/ee;
-			}
-		}
-	}
-
-	for (int i = 0; i < a_s; i++)
-	{
-		ee = 0.0;
-		for (int l = 0; l < len; l++)
-		{
-			ee += e[l + i*len] * e[l + i*len];
-		}
-		ee = sqrt(ee);
-
-		for (int l = 0; l < len; l++)
-		{
-			e[l + i*len] /= ee;
-		}
-	}
-	return;
-}
 
 double gctl::dist_inverse_weight(std::vector<double> *dis_vec, std::vector<double> *val_vec, int order)
 {

lib/algorithm/algorithm_func.h

@@ -0,0 +1,139 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_ALGORITHM_FUNC_H
+#define _GCTL_ALGORITHM_FUNC_H
+
+#include "../core/array.h"
+#include "../core/matrix.h"
+#include "../maths/mathfunc.h"
+
+namespace gctl
+{
+    /**
+	 * @brief      按距离反比加权计算均值
+	 *
+	 * @param      dis_vec  距离向量
+	 * @param      val_vec  数值向量
+	 * @param[in]  order    距离加权的阶次 默认为1
+	 *
+	 * @return     加权平均值
+	 */
+	double dist_inverse_weight(std::vector<double> *dis_vec, std::vector<double> *val_vec, int order = 1);
+
+	/**
+	 * @brief      查找一个数在已排序的数组中的位置，即找到包含该数的最小区间
+	 *
+	 * @param[in]  in_array    输入数组
+	 * @param[in]  array_size  数组大小
+	 * @param[in]  in_val      查找值
+	 * @param      index       返回的索引值
+	 *
+	 * @return     成功0失败-1
+	 */
+	int find_index(const double *in_array, int array_size, double in_val, int &index);
+
+	/**
+	 * @brief      查找一个数在已排序的数组中的位置，即找到包含该数的最小区间
+	 *
+	 * @param      in_array  输入数组
+	 * @param[in]  in_val    查找值
+	 * @param      index     返回的索引值
+	 *
+	 * @return     成功0失败-1
+	 */
+	int find_index(array<double> *in_array, double in_val, int &index);
+
+	/**
+	 * @brief      计算一维分形模型
+	 *
+	 * @param      out_arr     输出数组
+	 * @param[in]  l_val       分形计算的左端点值
+	 * @param[in]  r_val       分形计算的右端点值
+	 * @param[in]  maxi_range  最大变化值
+	 * @param[in]  smoothness  变化光滑度
+	 */
+	void fractal_model_1d(array<double> &out_arr, int out_size, double l_val, 
+		double r_val, double maxi_range, double smoothness);
+
+	/**
+	 * @brief      计算二维分形模型
+	 *
+	 * @param      out_arr     输出数组
+	 * @param[in]  dl_val      分形计算的左下角端点值
+	 * @param[in]  dr_val      分形计算的右下角端点值
+	 * @param[in]  ul_val      分形计算的左上角端点值
+	 * @param[in]  ur_val      分形计算的右上角端点值
+	 * @param[in]  maxi_range  最大变化值
+	 * @param[in]  smoothness  变化光滑度
+	 */
+	void fractal_model_2d(_2d_matrix &out_arr, int r_size, int c_size, double dl_val, 
+		double dr_val, double ul_val, double ur_val, double maxi_range, double smoothness, unsigned int seed = 0);
+
+    /**
+	 * @brief      一维数组差分（使用二阶差分公式）
+	 * 
+	 * 计算一个一维数组中相邻元素间的差分结果（求导）。
+	 *
+	 * @param[in]  in       输入数组
+	 * @param      diff     输出的差分结果
+	 * @param[in]  spacing  相邻元素的距离
+	 * @param[in]  order    求导的次数。最小为1（默认），最大为4。两边的数据将分别使用对应的向前或向后差分公式
+	 */
+	void difference_1d(const array<double> &in, array<double> &diff, double spacing, int order = 1);
+
+	/**
+	 * @brief      二维数组差分（使用二阶差分公式）
+	 * 
+	 * 计算一个二维数组中相邻元素间的差分结果（求导）。
+	 *
+	 * @param[in]  in       输入数组
+	 * @param      diff     输出的差分结果
+	 * @param[in]  spacing  相邻元素对应方向的距离
+	 * @param[in]  d_type   求导的类型
+	 * @param[in]  order    求导的次数。最小为1（默认），最大为4。边缘的数据将分别使用对应的向前或向后差分公式
+	 */
+	void difference_2d(const _2d_matrix &in, _2d_matrix &diff, double spacing, gradient_type_e d_type, int order = 1);
+
+	/**
+	 * @brief      二维数组差分（使用二阶差分公式）
+	 * 
+	 * 计算一个二维数组中相邻元素间的差分结果（求导）。数组以列优先的方式储存为一个一维数组
+	 *
+	 * @param[in]  in       输入数组
+	 * @param      diff     输出的差分结果
+	 * @param[in]  row_size 数组二维排列的行数
+	 * @param[in]  col_size 数组二维排列的列数
+	 * @param[in]  spacing  相邻元素对应方向的距离
+	 * @param[in]  d_type   求导的类型
+	 * @param[in]  order    求导的次数。最小为1（默认），最大为4。边缘的数据将分别使用对应的向前或向后差分公式
+	 */
+	void difference_2d(const array<double> &in, array<double> &diff, int row_size, int col_size, 
+		double spacing, gradient_type_e d_type, int order = 1);
+}
+
+#endif // _GCTL_ALGORITHM_FUNC_H

lib/algorithm/autodiff.h

@@ -56,6 +56,6 @@ namespace gctl
     private:
         double val_, der_;
     };
-};
+}
 
 #endif // _GCTL_AUTODIFF_H

lib/algorithm/boxsort2d.h

@@ -28,6 +28,7 @@
 #ifndef _BOXSORT2D_H
 #define _BOXSORT2D_H
 
+#include "../core.h"
 #include "heap_sort.h"
 
 namespace gctl
@@ -485,6 +486,6 @@ namespace gctl
         }
         return;
     }
-};
+}
 
 #endif // _BOXSORT2D_H

lib/algorithm/boxsort_sph.h

@@ -28,6 +28,7 @@
 #ifndef _BOXSORT_SPH_H
 #define _BOXSORT_SPH_H
 
+#include "../core.h"
 #include "heap_sort.h"
 
 namespace gctl

lib/algorithm/emd.h

@@ -29,6 +29,8 @@
 #define _GCTL_EMD_H
 
 // library head file
+#include "../gctl_config.h"
+#include "../core/array.h"
 #include "../core/matrix.h"
 
 #ifdef GCTL_EEMD
@@ -44,7 +46,7 @@ namespace gctl
     void ceemdan1d(const array<double> &in, matrix<double> &out, size_t M, 
         unsigned int ensemble_size, double noise_strength, unsigned int S_number, 
         unsigned int num_siftings, unsigned long int rng_seed);
-};
+}
 
 #endif // GCTL_EEMD
 

lib/algorithm/extrapolate.h

@@ -28,7 +28,9 @@
 #ifndef _GCTL_EXTRAPOLATE_H
 #define _GCTL_EXTRAPOLATE_H
 
-#include "gmath.h"
+#include "../core/array.h"
+#include "../core/matrix.h"
+#include "../maths/mathfunc_t.h"
 
 namespace gctl
 {
@@ -103,6 +105,6 @@ namespace gctl
 	 */
 	void zeros_extrapolate_2d(const array<double> &in_arr, array<double> &out_arr, int in_row, int in_col, 
 		int &out_row, int &out_col, int &ori_row, int &ori_col);
-};
+}
 
 #endif // _GCTL_EXTRAPOLATE_H

lib/algorithm/fft_filter.h

@@ -28,8 +28,9 @@
 #ifndef _FFT_FILTER_H
 #define _FFT_FILTER_H
 
+#include "../gctl_config.h"
 #include "../core/array.h"
-#include "fft.h"
+#include "../maths/fft.h"
 #include "extrapolate.h"
 
 #ifdef GCTL_FFTW3
@@ -93,7 +94,7 @@ namespace gctl
         array<double> fid_;
         array<std::complex<double> > freqs_;
     };
-};
+}
 
 #endif // GCTL_FFTW3
 

lib/algorithm/fir_filter.h

@@ -28,6 +28,7 @@
 #ifndef _RIF_FILTER_H
 #define _RIF_FILTER_H
 
+#include "../core.h"
 #include "windowfunc.h"
 
 namespace gctl
@@ -97,6 +98,6 @@ namespace gctl
 
         int taps_;       ///< Number of taps of the filter
     };
-};
+}
 
 #endif // _RIF_FILTER_H

lib/algorithm/glni.cpp

@@ -66,7 +66,7 @@ void glni::coefficients::initiate(size_t size)
 #ifdef GCTL_CHECK_SIZE
     if (size == 0 || size > 20)
     {
-        throw std::invalid_argument("Invalid initiating order. From glni::coefficients::initiate(...).");
+        throw invalid_argument("Invalid initiating order. From glni::coefficients::initiate(...).");
     }
 #endif
 
@@ -355,7 +355,7 @@ double glni::coefficients::node(size_t idx)
 #ifdef GCTL_CHECK_SIZE
     if (idx >= size_)
     {
-        throw std::out_of_range("Invalid index. From gctl::glni::coefficients::node(...)");
+        throw out_of_range("Invalid index. From gctl::glni::coefficients::node(...)");
     }
 #endif
 
@@ -367,7 +367,7 @@ double glni::coefficients::weight(size_t idx)
 #ifdef GCTL_CHECK_SIZE
     if (idx >= size_)
     {
-        throw std::out_of_range("Invalid index. From gctl::glni::coefficients::weight(...)");
+        throw out_of_range("Invalid index. From gctl::glni::coefficients::weight(...)");
     }
 #endif
 
@@ -404,7 +404,7 @@ void glni::line::initiate(size_t size)
 #ifdef GCTL_CHECK_SIZE
     if (size == 0 || size > 20)
     {
-        throw std::invalid_argument("Invalid initiating order. From glni::line::initiate(...)");
+        throw invalid_argument("Invalid initiating order. From glni::line::initiate(...)");
     }
 #endif
 
@@ -419,13 +419,13 @@ double glni::line::integral(double low, double hig, void *att)
 #ifdef GCTL_CHECK_BOUNDER
     if (hig <= low)
     {
-        throw std::length_error("Invalid integration range. From glni::line::integral(...)");
+        throw length_error("Invalid integration range. From glni::line::integral(...)");
     }
 #endif
 
     if (!initialized_)
     {
-        throw std::runtime_error("The object is not initialized. From gctl::glni::line::integral(...).");
+        throw runtime_error("The object is not initialized. From gctl::glni::line::integral(...).");
     }
 
     double scl = (hig - low)/2.0;
@@ -468,7 +468,7 @@ void glni::quadrilateral::initiate(size_t size)
 #ifdef GCTL_CHECK_SIZE
     if (size == 0 || size > 20)
     {
-        throw std::invalid_argument("Invalid initiating order. From glni::quadrilateral::initiate(...)");
+        throw invalid_argument("Invalid initiating order. From glni::quadrilateral::initiate(...)");
     }
 #endif
 
@@ -483,13 +483,13 @@ double glni::quadrilateral::integral(double xlow, double xhig, double ylow, doub
 #ifdef GCTL_CHECK_BOUNDER
     if (xhig <= xlow || yhig <= ylow)
     {
-        throw std::length_error("Invalid integration range. From glni::quadrilateral::integral(...)");
+        throw length_error("Invalid integration range. From glni::quadrilateral::integral(...)");
     }
 #endif
 
     if (!initialized_)
     {
-        throw std::runtime_error("Un-initialized. From glni::quadrilateral::integral(...)");
+        throw runtime_error("Un-initialized. From glni::quadrilateral::integral(...)");
     }
 
     double scl = (xhig-xlow)*(yhig-ylow)/4.0;
@@ -536,7 +536,7 @@ void glni::triangle::initiate(size_t size)
 #ifdef GCTL_CHECK_SIZE
     if (size == 0 || size > 20)
     {
-        throw std::invalid_argument("Invalid initiating order. From glni::triangle::initiate(...)");
+        throw invalid_argument("Invalid initiating order. From glni::triangle::initiate(...)");
     }
 #endif
 
@@ -552,13 +552,13 @@ double glni::triangle::integral(double x1, double x2, double x3, double y1,
 #ifdef GCTL_CHECK_BOUNDER
     if ((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1) <= GCTL_ZERO) // 这里没有设置为0是为了检测三个点是否共线
     {
-        throw std::length_error("Invalid integration range. From glni::triangle::integral(...)");
+        throw length_error("Invalid integration range. From glni::triangle::integral(...)");
     }
 #endif
 
     if (!initialized_)
     {
-        throw std::runtime_error("Un-initialized. From glni::triangle::integral(...)");
+        throw runtime_error("Un-initialized. From glni::triangle::integral(...)");
     }
 
     double ksi, eta;
@@ -609,7 +609,7 @@ void glni::triangle_c::initiate(size_t size)
 #ifdef GCTL_CHECK_SIZE
     if (size == 0 || size > 20)
     {
-        throw std::invalid_argument("Invalid initiating order. From glni::triangle_c::initiate(...)");
+        throw invalid_argument("Invalid initiating order. From glni::triangle_c::initiate(...)");
     }
 #endif
 
@@ -623,15 +623,15 @@ std::complex<double> glni::triangle_c::integral(double x1, double x2, double x3,
     double y2, double y3, void *att)
 {
 #ifndef GCTL_CHECK_BOUNDER
-    if ((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1) <= 1e-30) // 这里没有设置为0是为了检测三个点是否共线
+    if ((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1) <= GCTL_ZERO) // 这里没有设置为0是为了检测三个点是否共线
     {
-        throw std::length_error("Invalid integration range. From glni::triangle::integral(...)");
+        throw length_error("Invalid integration range. From glni::triangle::integral(...)");
     }
 #endif
 
     if (!initialized_)
     {
-        throw std::runtime_error("Un-initialized. From glni::triangle::integral(...)");
+        throw runtime_error("Un-initialized. From glni::triangle::integral(...)");
     }
 
 	double ksi, eta;
@@ -683,7 +683,7 @@ void glni::cube::initiate(size_t size)
 #ifdef GCTL_CHECK_SIZE
     if (size == 0 || size > 20)
     {
-        throw std::invalid_argument("Invalid initiating order. From glni::cube::initiate(...)");
+        throw invalid_argument("Invalid initiating order. From glni::cube::initiate(...)");
     }
 #endif
 
@@ -699,13 +699,13 @@ double glni::cube::integral(double x1, double x2, double y1, double y2,
 #ifdef GCTL_CHECK_BOUNDER
     if (x2 <= x1 || y2 <= y1 || z2 <= z1)
     {
-        throw std::length_error("Invalid integration range. From glni::cube::integral(...)");
+        throw length_error("Invalid integration range. From glni::cube::integral(...)");
     }
 #endif
 
     if (!initialized_)
     {
-        throw std::runtime_error("Un-initialized. From glni::cube::integral(...)");
+        throw runtime_error("Un-initialized. From glni::cube::integral(...)");
     }
 
     double gx, gy, gz, inte = 0.0;
@@ -758,7 +758,7 @@ void glni::tetrahedron::initiate(size_t size)
 #ifdef GCTL_CHECK_SIZE
     if (size == 0 || size > 20)
     {
-        throw std::invalid_argument("Invalid initiating order. From glni::tetrahedron::initiate(...)");
+        throw invalid_argument("Invalid initiating order. From glni::tetrahedron::initiate(...)");
     }
 #endif
 
@@ -775,14 +775,14 @@ double glni::tetrahedron::integral(double x1, double x2, double x3, double x4,
 {
     if (!initialized_)
     {
-        throw std::runtime_error("Un-initialized. From glni::tetrahedron::integral(...)");
+        throw runtime_error("Un-initialized. From glni::tetrahedron::integral(...)");
     }
 
     double ksi, eta, zta;
     double inte = 0.0;
     double scl = 6.0*tet_volume(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4);
 
-	if (scl < 1e-30)
+	if (scl < GCTL_ZERO)
     {
         throw std::invalid_argument("glni::tetrahedron::integral(...) got an invalid integration range.");
     }
@@ -835,7 +835,7 @@ void glni::triprism::initiate(size_t size)
 #ifdef GCTL_CHECK_SIZE
 	if (size == 0 || size > 20)
     {
-		throw std::std::invalid_argument("glni::triprism::initiate(...) got an invalid initiating order.");
+		throw std::invalid_argument("glni::triprism::initiate(...) got an invalid initiating order.");
     }
 #endif // GCTL_CHECK_BOUNDER
 

lib/algorithm/glni.h

@@ -28,7 +28,7 @@
 #ifndef _GCTL_GLNI_H
 #define _GCTL_GLNI_H
 
-#include <complex>
+#include "../core.h"
 #include "../gctl_config.h"
 
 namespace gctl
@@ -374,6 +374,6 @@ namespace gctl
             size_t size_;
         };
     }
-};
+}
 
 #endif // _GCTL_GLNI_H

lib/algorithm/heap_sort.h

@@ -29,7 +29,7 @@
 #define _GCTL_HEAPSORT_H
 
 // library's head file
-#include "../core/array.h"
+#include "../core.h"
 
 namespace gctl
 {
@@ -246,6 +246,6 @@ namespace gctl
 			return;
 		}
 	};
-};
+}
 
 #endif //_GCTL_HEAPSORT_H

lib/algorithm/interpolate.cpp

@@ -509,8 +509,8 @@ gctl::array<double> *gctl::p2p_dist_sph(point3ds * out_posi, int out_num, point3
 						{
 							tmp_id = box_list[search_id][b];
 							// 如果点到节点的距离小于椭球的半径则添加到向量中
-							tmp_pc1 = s2c(out_posi[i]);
+							tmp_pc1 = out_posi[i].s2c();
-							tmp_pc2 = s2c(in_posi[tmp_id]);
+							tmp_pc2 = in_posi[tmp_id].s2c();
 							tmp_dist = distance(tmp_pc1, tmp_pc2);
 
 							if (geometry3d::angle(tmp_pc1, tmp_pc2) <= search_deg && 

lib/algorithm/interpolate.h

@@ -28,7 +28,13 @@
 #ifndef _GCTL_INTERPOLATE_H
 #define _GCTL_INTERPOLATE_H
 
-#include "geometry3d.h"
+#include "../core/array.h"
+#include "../maths/mathfunc_t.h"
+#include "../maths/mathfunc.h"
+#include "../geometry/point2c.h"
+#include "../geometry/point3c.h"
+#include "../geometry/geometry3d.h"
+#include "algorithm_func.h"
 
 namespace gctl
 {
@@ -310,6 +316,6 @@ namespace gctl
 	 */
 	array<double> *p2p_dist_sph(point3ds * out_posi, int out_num, point3ds *in_posi, double *in_val, 
 		int in_num, double search_rlen, double search_deg);
-};
+}
 
 #endif // _GCTL_INTERPOLATE_H

lib/algorithm/kde.h

@@ -28,7 +28,9 @@
 #ifndef _GCTL_KERNEL_DENSITY_ESTIMATION_H
 #define _GCTL_KERNEL_DENSITY_ESTIMATION_H
 
-#include "../core/array.h"
+#include "../core.h"
+#include "../maths.h"
+#include "../io.h"
 
 namespace gctl
 {
@@ -169,6 +171,6 @@ namespace gctl
 
         double gaussian_kernel(double x, double y);
     };
-};
+}
 
 #endif // _GCTL_KERNEL_DENSITY_ESTIMATION_H

lib/algorithm/multinary.h

@@ -56,6 +56,6 @@ namespace gctl
         array<double> std_;
         array<double> ys_;
     };
-};
+}
 
 #endif // _GCTL_MULTINARY_H

lib/algorithm/sinkhorn.h

@@ -28,6 +28,8 @@
 #ifndef _GCTL_SINKHORN_H
 #define _GCTL_SINKHORN_H
 
+#include "../core/array.h"
+#include "algorithm_func.h"
 #include "interpolate.h"
 
 namespace gctl
@@ -112,6 +114,6 @@ namespace gctl
         matrix<double> RP_; // 整理后的转换矩阵
         matrix<double> rp_maxi_; // RP_中每一快的最大值
     };
-};
+}
 
 #endif // _GCTL_SINKHORN_H

lib/algorithm/space_filter.h

@@ -28,8 +28,8 @@
 #ifndef _GCTL_SPACE_FILTER_H
 #define _GCTL_SPACE_FILTER_H
 
-#include "../core/matrix.h"
+#include "../core.h"
-#include "mathfunc_ext.h"
+#include "../maths.h"
 
 namespace gctl
 {
@@ -161,6 +161,6 @@ namespace gctl
      * @param m_OrderNumClm 滑动窗口的列的多项式阶次
      */
     void trend_2d(matrix<double> &data, int m_RowNum, int m_ClmNum, int m_OrderNumRow, int m_OrderNumClm);
-};
+}
 
 #endif // _GCTL_SPACE_FILTER_H

lib/algorithm/variogram.h

@@ -137,6 +137,6 @@ namespace gctl
      * @return     Semi-variogram value
      */
     double variogram(const double &d, const variogram_model &vm, variogram_type_e type);
-};
+}
 
 #endif // _VARIOGRAM_H

lib/algorithm/windowfunc.h

@@ -28,7 +28,7 @@
 #ifndef _WINDOWFUNC_H
 #define _WINDOWFUNC_H
 
-#include "../core/array.h"
+#include "../core.h"
 
 namespace gctl
 {
@@ -63,6 +63,6 @@ namespace gctl
      * @param coeff Window coefficients
      */
     void window_blackman(int taps, array<double> &coeff);
-};
+}
 
 #endif // _WINDOWFUNC_H

lib/core.h

@@ -30,13 +30,12 @@
 
 #include "core/enum.h"
 #include "core/macro.h"
-#include "core/sptr.h"
 #include "core/vector_t.h"
 #include "core/exceptions.h"
+#include "core/sptr.h"
 #include "core/array.h"
 #include "core/matrix.h"
 #include "core/spmat.h"
 #include "core/sparray.h"
-#include "core/str.h"
 
 #endif // _GCTL_CORE_H

lib/core/array.h

@@ -28,15 +28,14 @@
 #ifndef _GCTL_ARRAY_H
 #define _GCTL_ARRAY_H
 
-#include <random>
-#include <chrono>
-
 // library's head files
 #include "../gctl_config.h"
 #include "enum.h"
 #include "macro.h"
 #include "vector_t.h"
 #include "exceptions.h"
+#include "random"
+#include "chrono"
 
 #ifdef GCTL_EIGEN
 /*The followings bypass a bug that VsCode's IntelliSense reports incorrect errors when using eigen3 library*/
@@ -63,8 +62,8 @@ namespace gctl
 	typedef array<float> _1f_array; ///< 1-D single-precision floating-point array
 	typedef array<double> _1d_array; ///< 1-D double-precision floating-point array
 	typedef array<std::string> _1s_array; ///< 1-D string array
-	typedef array<std::complex<float> > _1cf_array; ///< 1-D complex float array
+	typedef array<std::complex<float>> _1cf_array; ///< 1-D complex float array
-	typedef array<std::complex<double> > _1cd_array; ///< 1-D complex double array
+	typedef array<std::complex<double>> _1cd_array; ///< 1-D complex double array
 
 	/**
 	 * @brief Iterator template. This could be used to enable the use of iteration algorithms for array-like objects.
@@ -108,7 +107,7 @@ namespace gctl
 		 *
 		 * @param[in]  len   Length of the array. Must be equal to or bigger than zero.
 		 */
-		explicit array(size_t len);
+		array(size_t len);
 
 		/**
 		 * @brief      Construct an array with the given length and initial values.
@@ -146,7 +145,7 @@ namespace gctl
 		 * 
 		 * @param init_val Initial values
 		 */
-		explicit array(std::initializer_list<ArrValType> init_val);
+		array(std::initializer_list<ArrValType> init_val);
 
 		/**
 		 * @brief Construct a new array object as a sequence
@@ -584,6 +583,12 @@ namespace gctl
 		 */
 		void input(const Eigen::VectorXd &b);
 #endif // GCTL_EIGEN
+
+		/**
+		 * @brief element operate function pointer
+		 * 
+		 */
+		typedef void (*foreach_a_ptr)(ArrValType &ele_ptr, size_t id);
 		
 		/**
 		 * @brief 并行执行指定操作
@@ -774,7 +779,7 @@ namespace gctl
 		 * 
 		 * @return     返回的模长
 		 */
-		ArrValType module(norm_type_e n_type = L2) const;
+		ArrValType module(norm_type_e n_type = L2);
 
 		/**
 		 * @brief Normalize the array to the given module value
@@ -1504,6 +1509,18 @@ namespace gctl
 	}
 #endif // GCTL_EIGEN
 
+/*
+	template <typename ArrValType>
+	void array<ArrValType>::for_each(foreach_a_ptr func)
+	{
+		for (size_t i = 0; i < length_; i++)
+		{
+			func(val_[i], i);
+		}
+		return;
+	}
+*/
+
 	template <typename ArrValType>
 	void array<ArrValType>::show(std::ostream &os, char sep)
 	{
@@ -1623,7 +1640,7 @@ namespace gctl
 	{
 		if (rn*cn != length_)
 		{
-			throw std::invalid_argument("[gctl::array<T>::sequence2d] Invalid sequence sizes.");
+			throw invalid_argument("[gctl::array<T>::sequence2d] Invalid sequence sizes.");
 		}
 
 		array<ArrValType> out_x(cn), out_y(rn);
@@ -1647,7 +1664,7 @@ namespace gctl
 	{
 		if (ln*rn*cn != length_)
 		{
-			throw std::invalid_argument("[gctl::array<T>::sequence3d] Invalid sequence sizes.");
+			throw invalid_argument("[gctl::array<T>::sequence3d] Invalid sequence sizes.");
 		}
 
 		array<ArrValType> out_x(cn), out_y(rn), out_z(ln);
@@ -1692,7 +1709,7 @@ namespace gctl
 			throw std::runtime_error("Incompatible array sizes. gctl::array<T>::dot(...)");
 		}
 
-		ArrValType s = (ArrValType) 0;
+		ArrValType s = 0;
 		for (size_t i = 0; i < length_; i++)
 		{
 			s += val_[i]*a[i];
@@ -1727,32 +1744,15 @@ namespace gctl
 		return sum()/static_cast<ArrValType>(length_);
 	}
 
-	template <typename ArrValType>
-	ArrValType array<ArrValType>::variance() const
-	{
-		static_assert(std::is_arithmetic<ArrValType>::value, 
-			"gctl::array<T>::variance(...) could only be used with an arithmetic type.");
-
-		if (length_ == 0 || val_ == nullptr) return ArrValType{};
-
-		ArrValType m = mean();
-		ArrValType d = ArrValType(0);
-		for (size_t i = 0; i < length_; i++)
-		{
-			d += (val_[i] - m)*(val_[i] - m);
-		}
-		return d/static_cast<ArrValType>(length_);
-	}
-
 	template <typename ArrValType>
 	ArrValType array<ArrValType>::std() const
 	{
 		static_assert(std::is_arithmetic<ArrValType>::value, 
 			"gctl::array<T>::std(...) could only be used with an arithmetic type.");
 
-		if (length_ == 0 || val_ == nullptr) return ArrValType{};
+		if (length_ == 0 || length_ == 1) return ArrValType{};
 
-		ArrValType m = mean(), s = ArrValType(0);
+		ArrValType m = mean(), s = 0;
 		for (size_t i = 0; i < length_; i++)
 		{
 			s += (val_[i] - m)*(val_[i] - m);
@@ -1809,7 +1809,7 @@ namespace gctl
 	}
 
 	template <typename ArrValType>
-	ArrValType array<ArrValType>::module(norm_type_e n_type) const
+	ArrValType array<ArrValType>::module(norm_type_e n_type)
     {
         static_assert(std::is_arithmetic<ArrValType>::value, 
 			"gctl::array<T>::module(...) could only be used with an arithmetic type.");
@@ -1883,7 +1883,7 @@ namespace gctl
 
         if (length_ != b.size())
         {
-            throw std::runtime_error("[gctl::array<T>::orth] Incompatible array size.");
+            throw runtime_error("[gctl::array<T>::orth] Incompatible array size.");
         }
 
         ArrValType product = dot(b);
@@ -1913,13 +1913,28 @@ namespace gctl
         static_assert(std::is_arithmetic<ArrValType>::value, 
 			"gctl::array<T>::linear2log(...) could only be used with an arithmetic type.");
 
-        if (length_ == 0 || val_ == nullptr) return;
+        for (size_t i = 0; i < length_; ++i)
+        {
+            val_[i] = log(val_[i])/log(base);
+        }
+        return;
+    }
 
+    template <typename ArrValType>
+	ArrValType array<ArrValType>::variance() const
+    {
+        static_assert(std::is_arithmetic<ArrValType>::value, 
+			"gctl::array<T>::variance(...) could only be used with an arithmetic type.");
+
+		if (length_ == 0) return ArrValType{};
+
+        ArrValType mn = mean();
+        ArrValType d = 0;
         for (size_t i = 0; i < length_; i++)
         {
-            if (val_[i] <= ArrValType(0)) continue;  // 跳过非正值
+            d = d + (val_[i] - mn)*(val_[i] - mn);
-            val_[i] = std::log(val_[i])/std::log(base);
         }
+        return d/static_cast<ArrValType>(length_);
     }
 
     template <typename ArrValType>
@@ -1931,7 +1946,7 @@ namespace gctl
 		if (length_ < 2) return 0;
 
         ArrValType mn = mean();
-        ArrValType d = (ArrValType) 0;
+        ArrValType d = 0;
         for (size_t i = 0; i < length_; i++)
         {
             d = d + (val_[i] - mn)*(val_[i] - mn);
@@ -2029,6 +2044,6 @@ namespace gctl
 		else os << byte << " B\n";
 		return;
 	}
-};
+}
 
 #endif // _GCTL_ARRAY_H

lib/core/enum.h

@@ -257,6 +257,6 @@ namespace gctl
         SoftScale, // 所有数据按比例映射至max(min, min(arr))和min(max, max(arr))范围内
         CutOff, // 超过min和max范围数据直接设置为边界值
     };
-};
+}
 
 #endif // _GCTL_ENUM_H

lib/core/exceptions.h

@@ -28,11 +28,11 @@
 #ifndef _GCTL_EXCEPTIONS_H
 #define _GCTL_EXCEPTIONS_H
 
-#include <string>
+#include "string"
-#include <iostream>
+#include "iostream"
-#include <exception>
+#include "exception"
-#include <stdexcept>
+#include "stdexcept"
-#include <type_traits>
+#include "type_traits"
 
 /**
  * define error symbols
@@ -42,11 +42,11 @@
 #define GCTL_MESSAGE_ERROR -3
 
 #if defined _WINDOWS || __WIN32__
-#include <io.h>
+#include "io.h"
-#include <process.h>
+#include "process.h"
-#include <windows.h>
+#include "windows.h"
 #else
-#include <unistd.h>
+#include "unistd.h"
 #endif
 
 #if defined _WINDOWS || __WIN32__
@@ -113,28 +113,32 @@ namespace gctl
     {
     public:
         runtime_error() : std::runtime_error("GCTL: Unexpected runtime error."){}
-        explicit runtime_error(std::string estr) : std::runtime_error(("GCTL: Unexpected runtime error. "+estr).c_str()){}
+        runtime_error(std::string estr) : std::runtime_error(("GCTL: Unexpected runtime error. "+estr).c_str()){}
+        virtual ~runtime_error(){}
     };
 
     class range_error : public std::range_error
     {
     public:
         range_error() : std::range_error("GCTL: Invalid range detected."){}
-        explicit range_error(std::string estr) : std::range_error(("GCTL: Invalid range detected. "+estr).c_str()){}
+        range_error(std::string estr) : std::range_error(("GCTL: Invalid range detected. "+estr).c_str()){}
+        virtual ~range_error(){}
     };
 
     class overflow_error : public std::overflow_error
     {
     public:
         overflow_error() : std::overflow_error("GCTL: Overflow error."){}
-        explicit overflow_error(std::string estr) : std::overflow_error(("GCTL: Overflow error. "+estr).c_str()){}
+        overflow_error(std::string estr) : std::overflow_error(("GCTL: Overflow error. "+estr).c_str()){}
+        virtual ~overflow_error(){}
     };
 
     class underflow_error : public std::underflow_error
     {
     public:
         underflow_error() : std::underflow_error("GCTL: Underflow error."){}
-        explicit underflow_error(std::string estr) : std::underflow_error(("GCTL: Underflow error. "+estr).c_str()){}
+        underflow_error(std::string estr) : std::underflow_error(("GCTL: Underflow error. "+estr).c_str()){}
+        virtual ~underflow_error(){}
     };
 
 
@@ -142,36 +146,41 @@ namespace gctl
     {
     public:
         logic_error() : std::logic_error("GCTL: Logic error."){}
-        explicit logic_error(std::string estr) : std::logic_error(("GCTL: Logic error. "+estr).c_str()){}
+        logic_error(std::string estr) : std::logic_error(("GCTL: Logic error. "+estr).c_str()){}
+        virtual ~logic_error(){}
     };
 
     class domain_error : public std::domain_error
     {
     public:
         domain_error() : std::domain_error("GCTL: Domain error."){}
-        explicit domain_error(std::string estr) : std::domain_error(("GCTL: Domain error. "+estr).c_str()){}
+        domain_error(std::string estr) : std::domain_error(("GCTL: Domain error. "+estr).c_str()){}
+        virtual ~domain_error(){}
     };
 
     class invalid_argument : public std::invalid_argument
     {
     public:
         invalid_argument() : std::invalid_argument("GCTL: Invalid argument."){}
-        explicit invalid_argument(std::string estr) : std::invalid_argument(("GCTL: Invalid argument. "+estr).c_str()){}
+        invalid_argument(std::string estr) : std::invalid_argument(("GCTL: Invalid argument. "+estr).c_str()){}
+        virtual ~invalid_argument(){}
     };
 
     class length_error : public std::length_error
     {
     public:
         length_error() : std::length_error("GCTL: Length error."){}
-        explicit length_error(std::string estr) : std::length_error(("GCTL: Length error. "+estr).c_str()){}
+        length_error(std::string estr) : std::length_error(("GCTL: Length error. "+estr).c_str()){}
+        virtual ~length_error(){}
     };
 
     class out_of_range : public std::out_of_range
     {
     public:
         out_of_range() : std::out_of_range("GCTL: Out of range."){}
-        explicit out_of_range(std::string estr) : std::out_of_range(("GCTL: Out of range. "+estr).c_str()){}
+        out_of_range(std::string estr) : std::out_of_range(("GCTL: Out of range. "+estr).c_str()){}
+        virtual ~out_of_range(){}
     };
-};
+}
 
 #endif // _GCTL_EXCEPTIONS_H

lib/core/macro.h

@@ -43,13 +43,12 @@
 /**
  * physical constants
  */
-#define GCTL_WGS84_PoleRadius 6356752.3142 ///< WGS84椭球极半径
+#define GCTL_WGS84_PoleRadius 6356752.314 ///< WGS84椭球极半径
-#define GCTL_WGS84_EquatorRadius 6378136.460 ///< WGS84椭球长半径
+#define GCTL_WGS84_EquatorRadius 6378137 ///< WGS84椭球长半径
-#define GCTL_WGS84_FLAT 3.35281066474748e-3 ///<  (1.0/298.257223563)
+#define GCTL_WGS84_FLAT 3.35281066474748e-3 ///<  (1/298.25722356300)
 #define GCTL_WGS84_NormPotential 6.263685171456948e+07 ///< m**2/s**2
 #define GCTL_Earth_GM 3.986004418e+14 ///< m**3/s**2
 #define GCTL_Earth_Radius 6371008.8 ///< 地球平均半径
-#define GCTL_Earth_RefRadius 6371200.0 ///< 地球参考半径（常用于地磁场建模中）
 #define GCTL_Moon_Radius 1738000 ///< 月球平均半径
 // Ardalan A. A., Karimi R. & Grafarend E. W. (2010). A new reference equipotential surface, and reference ellipsoid for the planet mars.
 // Earth Moon Planet, 106:1-13

lib/core/matrix.h

@@ -28,22 +28,25 @@
 #ifndef _GCTL_MATRIX_H
 #define _GCTL_MATRIX_H
 
-#include <typeinfo>
 #include "array.h"
 
+#include "random"
+#include "chrono"
+#include "typeinfo"
+
 namespace gctl
 {
 	template <typename MatValType> class matrix;
 
 	/*
-	 * Here we define some commonly used matrix types.
+	 * Here we define some commonly used array types.
 	 */
-	typedef matrix<int> _2i_matrix; ///< 2-D integer matrix
+	typedef matrix<int> _2i_matrix; ///< 2-D integer array
-	typedef matrix<float> _2f_matrix; ///< 2-D single-precision floating-point matrix
+	typedef matrix<float> _2f_matrix; ///< 2-D single-precision floating-point array
-	typedef matrix<double> _2d_matrix; ///< 2-D double-precision floating-point matrix
+	typedef matrix<double> _2d_matrix; ///< 2-D double-precision floating-point array
-	typedef matrix<std::string> _2s_matrix; ///< 2-D string matrix
+	typedef matrix<std::string> _2s_matrix; ///< 2-D string array
-	typedef matrix<std::complex<float> > _2cf_matrix; ///< 1-D complex float matrix
+	typedef matrix<std::complex<float>> _2cf_matrix; ///< 1-D complex float array
-	typedef matrix<std::complex<double> > _2cd_matrix; ///< 1-D complex double matrix
+	typedef matrix<std::complex<double>> _2cd_matrix; ///< 1-D complex double array
 
 	/**
 	 * @brief      二维数组模版
@@ -114,7 +117,7 @@ namespace gctl
 		 *
 		 * @param[in]  b     输入向量的引用
 		 */
-		explicit matrix(const std::vector<std::vector<MatValType> > &b);
+		matrix(const std::vector<std::vector<MatValType> > &b);
 
 		/**
 		 * @brief      构造函数：从一维数组初始化二维数组
@@ -245,22 +248,6 @@ namespace gctl
 		 */
 		void assign_all(MatValType in_val);
 
-		/**
-		 * @brief 填充矩阵行
-		 *
-		 * @param r_id 行索引
-		 * @param in_val 行数据
-		 */
-		void fill_row(size_t r_id, const array<MatValType> &in_arr);
-
-		/**
-		 * @brief 填充矩阵列
-		 *
-		 * @param c_id 列索引
-		 * @param in_val 列数据
-		 */
-		void fill_column(size_t c_id, const array<MatValType> &in_arr);
-
 		/**
 		 * @brief      对全体元素进行缩放
 		 * 
@@ -275,20 +262,13 @@ namespace gctl
 		 * @param odr  计算行或列的平均值
 		 */
 		void mean(array<MatValType> &m, matrix_order_e odr = RowMajor);
-
-		/**
-		 * @brief      按行或列归一化矩阵元素
-		 *
-		 * @param odr  计算行或列的平均值
-		 */
-		void normalize(matrix_order_e odr = RowMajor);
 		
 		/**
 		 * @brief      计算矩阵的模长
 		 * 
 		 * @param n_type 模的类型
 		 */
-		double module(norm_type_e n_type = L2) const;
+		double norm(norm_type_e n_type = L2) const;
 
 		/**
 		 * @brief Set elements' value as a sequent.
@@ -806,32 +786,6 @@ namespace gctl
 		return;
 	}
 
-	template <typename MatValType>
-	void matrix<MatValType>::fill_row(size_t r_id, const array<MatValType> &in_arr)
-	{
-		if (r_id >= row_length) throw std::out_of_range("[gctl::matrix<T>::fill_row] Row index out of range.");
-		if (in_arr.size() != col_length) throw std::invalid_argument("[gctl::matrix<T>::fill_row] Invalid array length.");
-
-		for (size_t j = 0; j < col_length; j++)
-		{
-			val[r_id][j] = in_arr[j];
-		}
-		return;
-	}
-
-	template <typename MatValType>
-	void matrix<MatValType>::fill_column(size_t c_id, const array<MatValType> &in_arr)
-	{
-		if (c_id >= col_length) throw std::out_of_range("[gctl::matrix<T>::fill_column] Column index out of range.");
-		if (in_arr.size()!= row_length) throw std::invalid_argument("[gctl::matrix<T>::fill_column] Invalid array length.");
-
-		for (size_t i = 0; i < row_length; i++)
-		{
-			val[i][c_id] = in_arr[i];
-		}
-		return;
-	}
-
 	template <typename MatValType>
 	void matrix<MatValType>::scale(MatValType s)
 	{
@@ -878,47 +832,7 @@ namespace gctl
 	}
 
 	template <typename MatValType>
-	void matrix<MatValType>::normalize(matrix_order_e odr)
+	double matrix<MatValType>::norm(norm_type_e n_type) const
-	{
-		if (odr == RowMajor)
-		{
-			for (size_t i = 0; i < row_length; i++)
-			{
-				MatValType sum = 0.0;
-				for (size_t j = 0; j < col_length; j++)
-				{
-					sum += val[i][j]*val[i][j];
-				}
-				sum = sqrt(sum);
-
-				for (size_t j = 0; j < col_length; j++)
-				{
-					val[i][j] /= sum;
-				}
-			}
-		}
-		else // ColMajor
-		{
-			for (size_t j = 0; j < col_length; j++)
-			{
-				MatValType sum = 0.0;
-				for (size_t i = 0; i < row_length; i++)
-				{
-					sum += val[i][j]*val[i][j];
-				}
-				sum = sqrt(sum);
-
-				for (size_t i = 0; i < row_length; i++)
-				{
-					val[i][j] /= sum;
-				}
-			}
-		}
-		return;
-	}
-
-	template <typename MatValType>
-	double matrix<MatValType>::module(norm_type_e n_type) const
 	{
 		MatValType nval = (MatValType) 0;
 		if (n_type == L0)
@@ -1182,6 +1096,6 @@ namespace gctl
 		}
 		return;
 	}
-};
+}
 
 #endif // _GCTL_MATRIX_H

lib/core/sparray.h

@@ -30,6 +30,7 @@
 
 // library's head files
 #include "array.h"
+#include "enum.h"
 
 namespace gctl
 {
@@ -433,7 +434,7 @@ namespace gctl
 
 		for (int i = 0; i < arr1d.size(); i++)
 		{
-			b.set(arr1d[i].id, multipler*arr1d[i].val);
+			b.insert(arr1d[i].id, multipler*arr1d[i].val);
 		}
 		return;
 	}
@@ -653,6 +654,6 @@ namespace gctl
 	{
 		return arr2d.size();
 	}
-};
+}
 
 #endif // _GCTL_SPARRAY_H

lib/core/spmat.h

@@ -28,6 +28,8 @@
 #ifndef _GCTL_SPMAT_H
 #define _GCTL_SPMAT_H
 
+#include "enum.h"
+#include "array.h"
 #include "matrix.h"
 
 namespace gctl
@@ -40,8 +42,8 @@ namespace gctl
 	typedef spmat<int> _1i_spmat; ///< 整形浮点稀疏矩阵
 	typedef spmat<float> _1f_spmat; ///< 单精度浮点稀疏矩阵
 	typedef spmat<double> _1d_spmat; ///< 双精度浮点稀疏矩阵
-	typedef spmat<std::complex<float> > _1cf_spmat; ///< 单精度复浮点稀疏矩阵
+	typedef spmat<std::complex<float>> _1cf_spmat; ///< 单精度复浮点稀疏矩阵
-	typedef spmat<std::complex<double> > _1cd_spmat; ///< 双精度复浮点稀疏矩阵
+	typedef spmat<std::complex<double>> _1cd_spmat; ///< 双精度复浮点稀疏矩阵
 
 	/**
 	 * @brief      稀疏矩阵的节点元素类型
@@ -61,7 +63,6 @@ namespace gctl
 		 */
 		mat_node()
 		{
-			val = (NodeValueType) 0;
 			r_id = c_id = -1;
 			next_right = next_down = nullptr;
 		}
@@ -627,7 +628,6 @@ namespace gctl
 		r_end = c_end = nullptr;
 		d_link = nullptr;
 		r_count = c_count = nullptr;
-		zero_val = (MatrixValueType) 0;
 	}
 
 	template <typename MatrixValueType>
@@ -2941,7 +2941,6 @@ namespace gctl
 	{
 		if (!empty()) clear();
 		
-		zero_val = b.zero_val;
 		r_num = b.row_size();
 		c_num = b.col_size();
 		r_head = new mat_node<MatrixValueType>* [r_num];
@@ -3305,6 +3304,6 @@ namespace gctl
 
 #endif // GCTL_EIGEN
 
-};
+}
 
 #endif //_GCTL_SPMAT_H

lib/core/sptr.h

@@ -41,7 +41,7 @@ namespace gctl
 		//该类成员访问权限全部为private，因为不想让用户直接使用该类
 		//定义智能指针类为友元，因为智能指针类需要直接操纵辅助类
 		friend class smart_ptr<T>;
-		ref_ptr() : p(nullptr), count(0) {}
+		ref_ptr() : p(nullptr) {}
 		//构造函数的参数为基础对象的指针
 		ref_ptr(T *ptr) : p(ptr), count(1) {}
 		//析构函数
@@ -94,6 +94,6 @@ namespace gctl
 	private:
 		ref_ptr<T> *rp; //辅助类对象指针
 	};
-};
+}
 
 #endif //_GCTL_SPTR_H

lib/core/str.cpp

@@ -1,207 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#include "str.h"
-
-//替换str中所有lod_value为new_value,返回被替换的old_value的个数
-int gctl::replace_all(std::string& new_str, const std::string& old_str,const std::string& old_value,const std::string& new_value)
-{
-	int count = 0;
-	new_str = old_str;
-	for(std::string::size_type pos(0);pos!=std::string::npos;pos+=new_value.length())
-	{
-		if((pos=new_str.find(old_value,pos))!=std::string::npos)
-		{
-			new_str.replace(pos,old_value.length(),new_value);
-			count++;
-		}
-		else break;
-	}
-	return count;
-}
-
-//在输入字符串末尾添加一段字符串，如果输入字符串的尾端与待添加的字符串一致则不添加并返回原字符串
-std::string gctl::patch_string(std::string in_str, std::string patch_str)
-{
-	int inlen = in_str.length();
-	int exlen = patch_str.length();
-
-	std::string out_str = in_str;
-
-	if (exlen > inlen)
-	{
-		out_str += patch_str;
-		return out_str;
-	}
-
-	if (exlen == inlen && in_str != patch_str)
-	{
-		out_str += patch_str;
-		return out_str;
-	}
-	
-    if (in_str.substr(inlen - exlen, inlen) != patch_str)
-    {
-        out_str += patch_str;
-		return out_str;
-    }
-
-	return out_str;
-}
-
-//转换string对象为stringstream对象
-void gctl::str2ss(std::string in_str, std::stringstream &out_ss, std::string delimiter)
-{
-	if (delimiter != "")
-	{
-		std::string replace_str;
-		replace_all(replace_str, in_str, delimiter, " ");
-		out_ss.clear();
-		out_ss.str(replace_str);
-	}
-	else
-	{
-		out_ss.clear();
-		out_ss.str(in_str);
-	}
-	return;
-}
-
-/**
- * @brief      转换string字符串为double类型的数值
- * 
- * 这个函数的主要作用是检查输入字符串是否为nan或者inf等表示无效值的符号。有的编译器
- * 可以在>>输入符中完成此检测，但为了函数功能的稳定，所以在此处自定了这个函数。
- *
- * @param[in]  instr  输入字符串
- *
- * @return     返回的double类型的数值
- */
-double gctl::str2double(std::string instr)
-{
-	if (instr == "NAN" || instr == "nan" || instr == "NaN")
-		return NAN;
-	else if (instr == "INF" || instr == "inf" || instr == "Inf")
-		return INFINITY;
-	else
-	{
-		auto e(instr.find_first_of("dD"));
-		if (e != std::string::npos) instr[e] = 'e';
-
-		return atof(instr.c_str());
-	}
-}
-
-/**
- * @brief      转换double类型数值为string类型字符串
- *
- * @param[in]  in_d  输入数值
- *
- * @return     输出字符串
- */
-std::string gctl::double2str(double in_d)
-{
-	std::string tmp_str;
-	std::stringstream tmp_ss;
-
-	tmp_ss.clear();
-	if (std::isnan(in_d))
-		tmp_ss.str("nan");
-	else if (std::isinf(in_d))
-		tmp_ss.str("inf");
-	else tmp_ss << in_d;
-
-	tmp_ss >> tmp_str;
-	return tmp_str;
-}
-
-//返回指定长度的随机字符串
-void gctl::random_char(unsigned int length, char* out)
-{
-	char str[76] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz!@#$%^&*+_-=?";
-	int i,lstr;
-	char ss[2] = {0};
-	lstr = strlen(str);//计算字符串长度
-	srand((unsigned int)time((time_t *)NULL));//使用系统时间来初始化随机数发生器
-	//按指定大小返回相应的字符串
-	for(i = 1; i <= length; i++)
-	{
-		snprintf(ss, 2, "%c",str[(rand()%lstr)]);//rand()%lstr 可随机返回0-61之间的整数, str[0-61]可随机得到其中的字符
-		strcat(out,ss);//将随机生成的字符串连接到指定数组后面
-	}
-	return;
-}
-
-void gctl::random_str(unsigned int length, std::string &out_str)
-{
-	char *out_char = new char [length];
-	random_char(length, out_char);
-	out_str  = out_char;
-	delete[] out_char;
-	return;
-}
-
-void gctl::parse_string_with_quotes(std::string in_str, std::vector<std::string> &str_vec)
-{
-	std::vector<std::string> tmp_units;
-    parse_string_to_vector(in_str, ' ', tmp_units);
-
-    // Combine strings into one if they are enclosed by quotes
-    std::string tmp_str = "";
-    int i = 0;
-    while (i < tmp_units.size())
-    {
-        if (tmp_units[i].front() == '"' && tmp_units[i].back() == '"')
-        {
-            str_vec.push_back(tmp_units[i].substr(1, tmp_units[i].length() - 1));
-            i++;
-        }
-        else if (tmp_units[i].front() == '"')
-        {
-            tmp_str = tmp_units[i].substr(1, tmp_units[i].length());
-            i++;
-
-            while (tmp_units[i].back() != '"')
-            {
-                tmp_str += " " + tmp_units[i];
-                i++;
-            }
-
-            tmp_str += " " + tmp_units[i].substr(0, tmp_units[i].length() - 1);
-            i++;
-
-            str_vec.push_back(tmp_str);
-            tmp_str = "";
-        }
-        else
-        {
-            str_vec.push_back(tmp_units[i]);
-            i++;
-        }
-    }
-	return;
-}

lib/core/vector_t.h

@@ -28,9 +28,9 @@
 #ifndef _GCTL_VECTOR_TYPE_H
 #define _GCTL_VECTOR_TYPE_H
 
-#include <vector>
+#include "vector"
-#include <complex>
+#include "complex"
-#include <iostream>
+#include "iostream"
 
 namespace gctl
 {
@@ -41,14 +41,14 @@ namespace gctl
 	typedef std::vector<float> _1f_vector; ///< 1-D single-precision floating-point vector
 	typedef std::vector<double> _1d_vector; ///< 1-D double-precision floating-point vector
 	typedef std::vector<std::string> _1s_vector; ///< 1-D string vector
-	typedef std::vector<std::complex<float> > _1cf_vector; ///< 1-D complex float vector
+	typedef std::vector<std::complex<float>> _1cf_vector; ///< 1-D complex float vector
-	typedef std::vector<std::complex<double> > _1cd_vector; ///< 1-D complex double vector
+	typedef std::vector<std::complex<double>> _1cd_vector; ///< 1-D complex double vector
-	typedef std::vector<std::vector<int> > _2i_vector; ///< 2-D integer vector
+	typedef std::vector<std::vector<int>> _2i_vector; ///< 2-D integer vector
-	typedef std::vector<std::vector<float> > _2f_vector; ///< 2-D single-precision floating-point vector
+	typedef std::vector<std::vector<float>> _2f_vector; ///< 2-D single-precision floating-point vector
-	typedef std::vector<std::vector<double> > _2d_vector; ///< 2-D double-precision floating-point vector
+	typedef std::vector<std::vector<double>> _2d_vector; ///< 2-D double-precision floating-point vector
-	typedef std::vector<std::vector<std::string> > _2s_vector; ///< 2-D string vector
+	typedef std::vector<std::vector<std::string>> _2s_vector; ///< 2-D string vector
-	typedef std::vector<std::vector<std::complex<float> > > _2cf_vector; ///< 2-D complex float vector
+	typedef std::vector<std::vector<std::complex<float>>> _2cf_vector; ///< 2-D complex float vector
-	typedef std::vector<std::vector<std::complex<double> > > _2cd_vector; ///< 2-D complex double vector
+	typedef std::vector<std::vector<std::complex<double>>> _2cd_vector; ///< 2-D complex double vector
 
 	/**
 	 * @brief      Clear the memory space used by a 1-D vector
@@ -129,6 +129,6 @@ namespace gctl
 		while (os >> t){a.push_back(t);}
 		return os;
 	}
-};
+}
 
 #endif // _GCTL_VECTOR_TYPE_H

lib/geometry.h

@@ -1,4 +1,4 @@
-/********************************************************
+/********************************************************
  *  ██████╗  ██████╗████████╗██╗
  * ██╔════╝ ██╔════╝╚══██╔══╝██║
  * ██║  ███╗██║        ██║   ██║
@@ -25,11 +25,11 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#ifndef _GCTL_GRAPHIC_H
+#ifndef _GCTL_GEOMETRY_H
-#define _GCTL_GRAPHIC_H
+#define _GCTL_GEOMETRY_H
 
-#include "graphic/gnuplot.h"
+#include "geometry/refellipsoid.h"
-#include "graphic/gmt.h"
+#include "geometry/geometry2d.h"
-#include "graphic/cliplot.h"
+#include "geometry/geometry3d.h"
 
-#endif // _GCTL_GRAPHIC_H
+#endif // _GCTL_GEOMETRY_H

lib/geometry/block.h

@@ -274,6 +274,6 @@ namespace gctl
 		tar->ur = src->ur;
 		return;
 	}
-};
+}
 
 #endif // _GCTL_BLOCK_H

lib/geometry/edge.h

@@ -148,6 +148,6 @@ namespace gctl
         }
         return;
     }
-};
+}
 
 #endif // _GCTL_EDGE_H

lib/geometry/edge2d.h

@@ -28,6 +28,8 @@
 #ifndef _GCTL_EDGE2D_H
 #define _GCTL_EDGE2D_H
 
+#include "vertex.h"
+#include "entity.h"
 #include "triangle2d.h"
 
 namespace gctl
@@ -179,6 +181,6 @@ namespace gctl
         neigh[1] = &nei1;
         return;
     }
-};
+}
 
 #endif // _GCTL_EDGE2D_H

lib/geometry/entity.h

@@ -33,6 +33,7 @@
 // library head file
 #include "../core/array.h"
 
+
 namespace gctl
 {
     // this variable is only valid in this source file to control
@@ -250,6 +251,6 @@ namespace gctl
         }
         return;
     }
-};
+}
 
 #endif // _GCTL_ENTITY_H

lib/geometry/geometry2d.cpp

@@ -25,23 +25,21 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#include <cmath>
-#include <algorithm>
-#include <string>
-#include <iostream>
-#include <map>
-#include <vector>
-
 #include "geometry2d.h"
 
+#include "string"
+#include "cmath"
+#include "map"
+#include "algorithm"
+
 // 计算点t是否在直线ab上。
 bool gctl::geometry2d::collinear(const point2dc &a, const point2dc &b, 
 	const point2dc &t, double cut_off)
 {
 	point2dc at = t - a;
 	point2dc ab = b - a;
-	long double cos_arc = dot(at, ab)/(at.module()*ab.module());
+	double cos_arc = dot(at, ab)/(at.module()*ab.module());
-	if (sqrt(1.0 - cos_arc*cos_arc)*at.module() <= cut_off)
+	if (sqrt(1 - cos_arc*cos_arc)*at.module() <= cut_off)
 		return true;
 	else return false;
 }

lib/geometry/geometry2d.h

@@ -28,8 +28,14 @@
 #ifndef _GCTL_GEOMETRY2D_H
 #define _GCTL_GEOMETRY2D_H
 
-#include "../poly/triangle2d.h"
+// library's head files
-#include "../poly/edge2d.h"
+#include "../core/array.h"
+#include "../maths.h"
+
+#include "triangle2d.h"
+#include "triangle2d2o.h"
+#include "edge2d.h"
+#include "rectangle2d.h"
 
 namespace gctl
 {
@@ -233,6 +239,5 @@ namespace gctl
 		void extract_triangular_mesh_2d(const array<triangle2d> &in_eles, const array<bool> &out_list, 
 			array<vertex2dc> &out_nodes, array<triangle2d> &out_eles);
 	}
-};
+}
-
 #endif // _GCTL_GEOMETRY2D_H

lib/geometry/geometry3d.cpp

@@ -25,16 +25,15 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#include <algorithm>
-#include <string>
-#include <cmath>
-#include <vector>
-#include <map>
-#include <sstream>
-#include <iomanip>
-
 #include "geometry3d.h"
 
+#include "string"
+#include "cmath"
+#include "vector"
+#include "map"
+#include "algorithm"
+
+
 double gctl::geometry3d::angle(const point3dc &a, const point3dc &b) //向量夹角
 {
 	double d = GCTL_SetToBox(-1.0, 1.0, dot(a, b)/(a.module() * b.module()));
@@ -111,7 +110,7 @@ void gctl::geometry3d::get_plane_coeff(double x1, double x2, double x3, double y
 
 void gctl::geometry3d::get_plane_coeff(const point3ds &p1, const point3ds &p2, double &A, double &B, double &C, double &D)
 {
-	gctl::point3dc c1 = s2c(p1), c2 = s2c(p2);
+	gctl::point3dc c1 = p1.s2c(), c2 = p2.s2c();
 	get_plane_coeff(c1.x, c2.x, 0.0, c1.y, c2.y, 0.0, c1.z, c2.z, 0.0, A, B, C, D);
 	return;
 }
@@ -136,11 +135,11 @@ gctl::point3ds gctl::geometry3d::track_sphere_arc(const point3ds &v1, const poin
 		throw std::runtime_error("[gctl::geometry3d::track_sphere_arc] Invalid radius.");
 	}
 
-	point3dc p1 = s2c(v1), p2 = s2c(v2);
+	point3dc p1 = v1.s2c(), p2 = v2.s2c();
 	double Phi = angle(p1, p2); // Phi为弧度表示
 	point3dc pc = cos(arc)*p1 + sin(arc)/sin(Phi)*(p2 - cos(Phi)*p1);
 
-	return c2s(pc);
+	return pc.c2s();
 }
 
 gctl::point3dc gctl::geometry3d::dot_on_plane(const point3dc &face_dot, const point3dc &face_nor, 

lib/geometry/geometry3d.h

@@ -28,9 +28,21 @@
 #ifndef _GCTL_GEOMETRY3D_H
 #define _GCTL_GEOMETRY3D_H
 
-#include "../poly/triangle.h"
+// library's head file
-#include "../poly/tetrahedron.h"
+#include "../core/array.h"
-#include "gmath.h"
+#include "../core/matrix.h"
+#include "../maths.h"
+
+#include "tensor.h"
+#include "node.h"
+#include "edge.h"
+#include "tetrahedron.h"
+#include "triangle.h"
+#include "block.h"
+#include "tri_cone.h"
+#include "sphere.h"
+#include "tesseroid.h"
+#include "prism.h"
 
 namespace gctl
 {
@@ -349,6 +361,5 @@ namespace gctl
 		void cut_triangular_mesh(const array<triangle> &in_eles, const point3dc &nor, const point3dc &surf, 
 			array<vertex3dc> &out_nodes, array<triangle> &out_eles);
 	}
-};
+}
-
 #endif // _GCTL_GEOMETRY3D_H

lib/geometry/node.h

@@ -147,6 +147,6 @@ namespace gctl
 		copy_entity(tar, src);
 		return;
 	}
-};
+}
 
 #endif // _GCTL_NODE_H

lib/geometry/point2c.h

@@ -28,23 +28,30 @@
 #ifndef _GCTL_POINT2C_H
 #define _GCTL_POINT2C_H
 
+#include "../core/macro.h"
+#include "../core/exceptions.h"
+#include "point2p.h"
+
 #include "iostream"
 #include "string"
 #include "cmath"
 #include "iomanip"
 #include "regex"
 
-#include "../core/array.h"
-
 namespace gctl
 {
-	// static variable for controlling the IO process
+	template <typename T> struct point2c;
-	static int p2c_io_psn = 6;
+	template <typename T> struct point2p;
 
-	template <typename T> class point2c;
 	typedef point2c<double> point2dc;
 	typedef point2c<float>  point2fc;
 
+#ifndef IO_PSN
+#define IO_PSN
+	// static variable for controlling the IO process
+	static int io_psn = 6;
+#endif // IO_PSN
+
 	/**
 	 * @brief      A point under the 2D Cartesian coordinates (aka. 2D vector).
 	 */
@@ -118,7 +125,7 @@ namespace gctl
 		 * 
 		 * @return The point2dp object.
 		 */
-		//point2p<T> c2p() const;
+		point2p<T> c2p() const;
 		/**
 		 * @brief Set a point2c object from a string. The accepted format: (x, y)
 		 * 
@@ -246,7 +253,7 @@ namespace gctl
 	{
 		return point2c<T>(x/module(), y/module());
 	}
-/*
+
 	template <typename T>
 	gctl::point2p<T> gctl::point2c<T>::c2p() const
 	{
@@ -257,7 +264,7 @@ namespace gctl
 		else outp.arc = atan2(y, x) + 2.0*GCTL_Pi;
 		return outp;
 	}
-*/
+
 	template <typename T>
 	void gctl::point2c<T>::str(std::string str)
 	{
@@ -282,7 +289,7 @@ namespace gctl
 			throw invalid_argument("Invalid precision. From point2c::set_io_precision(...)");
 		}
 
-		p2c_io_psn = psn;
+		io_psn = psn;
 		return;
 	}
 
@@ -298,7 +305,7 @@ namespace gctl
 	template <typename T>
 	void gctl::point2c<T>::out_loc(std::ostream &os, char deli) const
 	{
-		os << std::setprecision(p2c_io_psn) << x << deli << y;
+		os << std::setprecision(io_psn) << x << deli << y;
 		return;
 	}
 
@@ -412,7 +419,7 @@ namespace gctl
 	template <typename T>
 	std::ostream &operator <<(std::ostream & os, const point2c<T> &a)
 	{
-		os << std::setprecision(p2c_io_psn) << a.x << " " << a.y;
+		os << std::setprecision(io_psn) << a.x << " " << a.y;
 		return os;
 	}
 
@@ -542,6 +549,6 @@ namespace gctl
 		}
 		return;
 	}
-};
+}
 
 #endif // _GCTL_POINT2C_H

lib/geometry/point2p.h

@@ -28,24 +28,30 @@
 #ifndef _GCTL_POINT2P_H
 #define _GCTL_POINT2P_H
 
+#include "../core/macro.h"
+#include "../core/exceptions.h"
+#include "point2c.h"
+
 #include "iostream"
 #include "string"
 #include "cmath"
 #include "iomanip"
 #include "regex"
 
-#include "../core/macro.h"
-#include "../core/exceptions.h"
-
 namespace gctl
 {
-	// static variable for controlling the IO process
+	template <typename T> struct point2c;
-	static int p2p_io_psn = 6;
+	template <typename T> struct point2p;
 
-	template <typename T> class point2p;
 	typedef point2p<double> point2dp;
 	typedef point2p<float>  point2fp;
 
+#ifndef IO_PSN
+#define IO_PSN
+	// static variable for controlling the IO process
+	static int io_psn = 6;
+#endif // IO_PSN
+
 	/**
 	 * @brief      A point under the 2-D Polar coordinates.
 	 */
@@ -110,7 +116,7 @@ namespace gctl
 		 * 
 		 * @return The point2c object.
 		 */
-		//point2c<T> p2c() const;
+		point2c<T> p2c() const;
 		/**
 		 * @brief      Parse the point's parameters from a string
 		 *
@@ -218,7 +224,7 @@ namespace gctl
 		else if (arc < 0.0) arc = fmod(arc, (2.0*GCTL_Pi)) + 2.0*GCTL_Pi;
 		return;
 	}
-	/*
+	
 	template <typename T>
 	gctl::point2c<T> gctl::point2p<T>::p2c() const
 	{
@@ -227,7 +233,7 @@ namespace gctl
 		outc.y = rad * sin(arc);
 		return outc;
 	}
-	*/
+
 	template <typename T>
 	void gctl::point2p<T>::str(std::string str)
 	{
@@ -253,14 +259,14 @@ namespace gctl
 			throw invalid_argument("Invalid precision. From point2p::set_io_precision(...)");
 		}
 
-		p2p_io_psn = psn;
+		io_psn = psn;
 		return;
 	}
 
 	template <typename T>
 	void gctl::point2p<T>::out_loc(std::ostream &os, char deli) const
 	{
-		os << std::setprecision(p2p_io_psn) << rad << deli << arc;
+		os << std::setprecision(io_psn) << rad << deli << arc;
 		return;
 	}
 	
@@ -302,7 +308,7 @@ namespace gctl
 	template <typename T>
 	std::ostream &operator <<(std::ostream & os, const point2p<T> &a)
 	{
-		os << std::setprecision(p2p_io_psn) << a.rad << " " << a.arc;
+		os << std::setprecision(io_psn) << a.rad << " " << a.arc;
 		return os;
 	}
 
@@ -383,6 +389,6 @@ namespace gctl
 		}
 		return;
 	}
-};
+}
 
 #endif // _GCTL_POINT2P_H

lib/geometry/point3c.h

@@ -28,23 +28,31 @@
 #ifndef _GCTL_POINT3C_H
 #define _GCTL_POINT3C_H
 
+#include "../core/macro.h"
+#include "../core/exceptions.h"
+#include "../core/array.h"
+#include "point3s.h"
+
 #include "iostream"
 #include "string"
 #include "cmath"
 #include "iomanip"
 #include "regex"
 
-#include "../core/array.h"
-
 namespace gctl
 {
-	// static variable for controlling the IO process
+	template <typename T> struct point3c;
-	static int p3c_io_psn = 6;
+	template <typename T> struct point3s;
 
-	template <typename T> class point3c;
 	typedef point3c<double> point3dc;
 	typedef point3c<float>  point3fc;
 
+#ifndef IO_PSN
+#define IO_PSN
+	// static variable for controlling the IO process
+	static int io_psn = 6;
+#endif // IO_PSN
+
 	/**
 	 * @brief      直角坐标系内的实点
 	 */
@@ -112,7 +120,7 @@ namespace gctl
 		 *
 		 * @return     spherical coordinates point
 		 */
-		//point3s<T> c2s() const;
+		point3s<T> c2s() const;
 		/**
 		 * @brief      Return the normal vector
 		 *
@@ -263,7 +271,7 @@ namespace gctl
 	{
 		return sqrt(x*x + y*y + z*z);
 	}
-/*
+
 	template <typename T>
 	gctl::point3s<T> gctl::point3c<T>::c2s() const
 	{
@@ -278,7 +286,7 @@ namespace gctl
 		outs.lon = atan2(y, x)*180.0/GCTL_Pi;
 		return outs;
 	}
-*/
+
 	template <typename T>
 	gctl::point3c<T> gctl::point3c<T>::normal() const
 	{
@@ -320,7 +328,7 @@ namespace gctl
 			throw invalid_argument("Invalid precision. From point3c::set_io_precision(...)");
 		}
 
-		p3c_io_psn = psn;
+		io_psn = psn;
 		return;
 	}
 
@@ -395,7 +403,7 @@ namespace gctl
 	template <typename T>
 	void gctl::point3c<T>::out_loc(std::ostream &os, char deli) const
 	{
-		os << std::setprecision(p3c_io_psn) << x << deli << y << deli << z;
+		os << std::setprecision(io_psn) << x << deli << y << deli << z;
 		return;
 	}
 
@@ -519,7 +527,7 @@ namespace gctl
 	template <typename T>
 	std::ostream &operator <<(std::ostream & os, const point3c<T> &a)
 	{
-		os << std::setprecision(p3c_io_psn) << a.x << " " << a.y << " " << a.z;
+		os << std::setprecision(io_psn) << a.x << " " << a.y << " " << a.z;
 		return os;
 	}
 
@@ -706,6 +714,6 @@ namespace gctl
 		}
 		return;
 	}
-};
+}
 
 #endif // _GCTL_POINT3C_H

lib/geometry/point3s.h

@@ -28,24 +28,30 @@
 #ifndef _GCTL_POINT3S_H
 #define _GCTL_POINT3S_H
 
+#include "../core/macro.h"
+#include "../core/exceptions.h"
+#include "point3c.h"
+
 #include "iostream"
 #include "string"
 #include "cmath"
 #include "iomanip"
 #include "regex"
 
-#include "../core/macro.h"
-#include "../core/exceptions.h"
-
 namespace gctl
 {
-	// static variable for controlling the IO process
+	template <typename T> struct point3c;
-	static int p3s_io_psn = 6;
+	template <typename T> struct point3s;
 
-	template <typename T> class point3s;
 	typedef point3s<double> point3ds;
 	typedef point3s<float>  point3fs;
 
+#ifndef IO_PSN
+#define IO_PSN
+	// static variable for controlling the IO process
+	static int io_psn = 6;
+#endif // IO_PSN
+
 	/**
 	 * @brief      球坐标系中的一个点（经纬度表示的）
 	 */
@@ -64,7 +70,7 @@ namespace gctl
 		void set(T r, T ln, T lt);
 		void set(const point3s<T> &b);
 		double module() const;
-		//point3c<T> s2c() const; // 转换为球坐标点
+		point3c<T> s2c() const; // 转换为球坐标点
 		void str(std::string str);
 		/**
 		 * @brief      Sets the i/o precision of the type.
@@ -172,15 +178,17 @@ namespace gctl
 	{
 		return rad;
 	}
-/*
+
 	template <typename T>
 	gctl::point3c<T> gctl::point3s<T>::s2c() const
 	{
-		//point3c<T> outc;
+		/*
-		//outc.x = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*cos((2.0 + lon/180.0)*GCTL_Pi);
+		point3c<T> outc;
-		//outc.y = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*sin((2.0 + lon/180.0)*GCTL_Pi);
+		outc.x = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*cos((2.0 + lon/180.0)*GCTL_Pi);
-		//outc.z = rad*cos((0.5 - lat/180.0)*GCTL_Pi);
+		outc.y = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*sin((2.0 + lon/180.0)*GCTL_Pi);
-		//return outc;
+		outc.z = rad*cos((0.5 - lat/180.0)*GCTL_Pi);
+		return outc;
+		*/
 
 		point3c<T> v;
 		v.x = rad*cos(GCTL_Pi*lat/180.0)*cos(GCTL_Pi*lon/180.0);
@@ -188,7 +196,7 @@ namespace gctl
 		v.z = rad*sin(GCTL_Pi*lat/180.0);
 		return v;
 	}
-*/
+
 	template <typename T>
 	void gctl::point3s<T>::str(std::string str)
 	{
@@ -214,14 +222,14 @@ namespace gctl
 			throw invalid_argument("Invalid precision. From point3s::set_io_precision(...)");
 		}
 
-		p3s_io_psn = psn;
+		io_psn = psn;
 		return;
 	}
 
 	template <typename T>
 	void gctl::point3s<T>::out_loc(std::ostream &os, char deli) const
 	{
-		os << std::setprecision(p3s_io_psn) << rad << deli << lon << deli << lat;
+		os << std::setprecision(io_psn) << rad << deli << lon << deli << lat;
 		return;
 	}
 
@@ -398,7 +406,7 @@ namespace gctl
 	template <typename T>
 	std::ostream &operator <<(std::ostream & os, const point3s<T> &a)
 	{
-		os << std::setprecision(p3s_io_psn) << a.rad << " " << a.lon << " " << a.lat;
+		os << std::setprecision(io_psn) << a.rad << " " << a.lon << " " << a.lat;
 		return os;
 	}
 
@@ -506,6 +514,6 @@ namespace gctl
 		}
 		return;
 	}
-};
+}
 
 #endif // _GCTL_POINT3S_H

lib/geometry/prism.h

@@ -178,6 +178,6 @@ namespace gctl
 		copy_entity(tar, src);
 		return;
 	}
-};
+}
 
 #endif // _GCTL_PRISM_H

lib/geometry/rectangle2d.h

@@ -204,6 +204,6 @@ namespace gctl
         }
         return;
     }
-};
+}
 
 #endif // _GCTL_RECTANGLE2D_H

lib/geometry/refellipsoid.cpp

@@ -25,43 +25,40 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#ifndef _GCTL_MATH_H
+#include "refellipsoid.h"
-#define _GCTL_MATH_H
 
-#include "math/gmath.h"
+gctl::refellipsoid::refellipsoid(/* args */){}
 
-#include "math/legendre.h"
+gctl::refellipsoid::refellipsoid(refellipsoid_type_e refellipsoid)
-#include "math/linear_algebra.h"
+{
+    set(refellipsoid);
+}
 
-#include "math/heap_sort.h"
+gctl::refellipsoid::refellipsoid(double R, double r, bool flat)
-#include "math/boxsort2d.h"
+{
-#include "math/boxsort_sph.h"
+    set(R, r, flat);
-#include "math/variogram.h"
+}
 
-#include "math/interpolate.h"
+gctl::refellipsoid::~refellipsoid(){}
-#include "math/extrapolate.h"
 
-#include "math/gaussfunc.h"
+void gctl::refellipsoid::set(refellipsoid_type_e refellipsoid)
-#include "math/sinkhorn.h"
+{
-#include "math/kde.h"
+   if (refellipsoid == Earth) {r_ = R_ = GCTL_Earth_Radius;}
-#include "math/autodiff.h"
+   else if (refellipsoid == Moon) {r_ = R_ = GCTL_Moon_Radius;}
-#include "math/multinary.h"
+   else if (refellipsoid == Mars) {r_ = R_ = GCTL_Mars_Radius;}
-#include "math/emd.h"
+   else if (refellipsoid == WGS84) {r_ = GCTL_WGS84_PoleRadius; R_ = GCTL_WGS84_EquatorRadius;}
+   else if (refellipsoid == Ardalan2010) {r_ = GCTL_Mars_Ardalan2010_PoleRadius; R_ = GCTL_Mars_Ardalan2010_EquatorRadius;}
+   else throw std::invalid_argument("Invalid reference system type for gctl::refellipsoid::set(...)");
+}
 
-#include "math/geometry2d.h"
+void gctl::refellipsoid::set(double R, double r_or_flat, bool is_flat)
-#include "math/geometry3d.h"
+{
-#include "math/refellipsoid.h"
+    R_ = R;
+    if (is_flat) r_ = R_*(1 - 1/r_or_flat);
+    else r_ = r_or_flat;
+}
 
-#include "math/mathfunc_ext.h"
+double gctl::refellipsoid::radius_at(double lati)
-#include "math/space_filter.h"
+{
-
+    return ellipse_radius_2d(R_, r_, lati*M_PI/180.0, 0.0);
-#include "math/windowfunc.h"
+}
-#include "math/fir_filter.h"
-
-#include "math/fft.h"
-#include "math/fft_filter.h"
-
-#include "math/shapefunc.h"
-#include "math/glni.h"
-
-#endif // _GCTL_MATH_H

lib/geometry/refellipsoid.h

@@ -28,8 +28,9 @@
 #ifndef _GCTL_REFELLIPSOID_H
 #define _GCTL_REFELLIPSOID_H
 
-#include "../poly/vertex.h"
+#include "../core/macro.h"
-#include "gmath.h"
+#include "../core/exceptions.h"
+#include "../maths/mathfunc.h"
 
 namespace gctl
 {
@@ -37,7 +38,6 @@ namespace gctl
     {
         // reference system named as the planet name represents a mean radius reference sphere
         Earth,
-        MagEarth,
         Moon,
         Mars,
         // Reference ellipsoids 
@@ -75,48 +75,17 @@ namespace gctl
         void set(double R, double r_or_flat, bool is_flat = false);
 
         /**
-         * @brief Compute the local rRdius of curvature on the reference ellipsoid
+         * @brief Get the ellipsoidal radius at the inquiring latitude 
          * 
-         * @param geodetic_lati the inquiring latitude
+         * @param lati the inquiring latitude
-         * @return radius
+         * @return ellipsoidal radius
          */
-        double geodetic_radius(double geodetic_lati);
+        double radius_at(double lati);
-
-        /**
-         * @brief Converts Geodetic coordinates to Spherical coordinates
-         * 
-         * @param geodetic_lati Geodetic latitude
-         * @param geodetic_hei Height above the ellipsoid
-         * @param sph_lati Spherical latitude
-         * @param sph_rad Spherical radius
-         */
-        void geodetic2spherical(double geodetic_lati, double geodetic_hei, 
-            double& sph_lati, double& sph_rad);
-
-        void spherical2geodetic(const point3ds& ps, double& geodetic_lon, double& geodetic_lati, 
-            double& geodetic_hei, double eps = 1e-16, int cnt = 30);
-
-        /**
-         * @brief Convert xyz coordinates to Geodetic coodinates
-         * 
-         * @param pc A point
-         * @param geodetic_lon Geodetic longitude
-         * @param geodetic_lati Geodetic latitude
-         * @param geodetic_hei Height above the ellipsoid
-         * @param eps solving precision
-         * @param cnt maximal iteration
-         */
-        void xyz2geodetic(const point3dc& pc, double& geodetic_lon, double& geodetic_lati, 
-            double& geodetic_hei, double eps = 1e-16, int cnt = 30);
-
-        void geodetic2xyz(double geodetic_lon, double geodetic_lati, 
-            double geodetic_hei, point3dc& pc);
 
     private:
-        double r_, R_, f_; // semi-minor, semi-major and flat rate
+        double r_, R_;
-        double e_;
     };
     
-};
+}
 
 #endif // _GCTL_REFELLIPSOID_H

lib/geometry/sphere.h

@@ -171,6 +171,6 @@ namespace gctl
         tar->rad = src->rad;
 		return;
 	}
-};
+}
 
 #endif // _GCTL_SPHERE_H

lib/geometry/tensor.h

@@ -28,7 +28,8 @@
 #ifndef _GCTL_TENSOR_H
 #define _GCTL_TENSOR_H
 
-#include "../core/matrix.h"
+#include "../core.h"
+
 #include "point3c.h"
 
 namespace gctl
@@ -396,6 +397,6 @@ namespace gctl
 		t.val[2][0]=a.z*b.x; t.val[2][1]=a.z*b.y; t.val[2][2]=a.z*b.z;
 		return t;
 	}
-};
+}
 
 #endif // _GCTL_TENSOR_H

lib/geometry/tesseroid.h

@@ -273,6 +273,6 @@ namespace gctl
 		tar->ur = src->ur;
 		return;
 	}
-};
+}
 
 #endif // _GCTL_TESSEROID_H

lib/geometry/tetrahedron.h

@@ -336,10 +336,10 @@ namespace gctl
         }
         this->self_host = true;
 
-        point3dc p0 = s2c(ps0);
+        point3dc p0 = ps0.s2c();
-        point3dc p1 = s2c(ps1);
+        point3dc p1 = ps1.s2c();
-        point3dc p2 = s2c(ps2);
+        point3dc p2 = ps2.s2c();
-        point3dc p3 = s2c(ps3);
+        point3dc p3 = ps3.s2c();
 
         this->id = index;
         this->vert[0]->set(p0, 4*index + 0);
@@ -577,6 +577,6 @@ namespace gctl
 		tar->vec_order = src->vec_order;
 		return;
 	}
-};
+}
 
 #endif // _GCTL_TETRAHEDRON_H

lib/geometry/tri_cone.h

@@ -179,9 +179,9 @@ namespace gctl
         }
         this->self_host = true;
 
-        point3dc p0 = s2c(ps0);
+        point3dc p0 = ps0.s2c();
-        point3dc p1 = s2c(ps1);
+        point3dc p1 = ps1.s2c();
-        point3dc p2 = s2c(ps2);
+        point3dc p2 = ps2.s2c();
 
         // 重新检查顶点的排序 初始化为逆时针排序
         point3dc tmp_p;

lib/geometry/triangle.h

@@ -221,6 +221,6 @@ namespace gctl
 		hgt_x = sqrt(edge_len[(c_id+2)%3]*edge_len[(c_id+2)%3] - hgt_y*hgt_y);
 		return;
 	}
-};
+}
 
 #endif // _GCTL_TRIANGLE_H

lib/geometry/triangle2d.h

@@ -214,6 +214,6 @@ namespace gctl
 		copy_entity(tar, src);
 		return;
 	}
-};
+}
 
 #endif // _GCTL_TRIANGLE2D_H

lib/geometry/triangle2d2o.h

@@ -296,6 +296,6 @@ namespace gctl
 		copy_entity(tar, src);
 		return;
 	}
-};
+}
 
 #endif // _GCTL_TRIANGLE2D2O_H

lib/geometry/vertex.h

@@ -35,57 +35,6 @@
 
 namespace gctl
 {
-    template <typename T>
-	gctl::point2p<T> c2p(const gctl::point2c<T>& c)
-	{
-		point2p<T> outp;
-		outp.rad = c.module();
-		if (c.y >= 0.0)
-			outp.arc= atan2(c.y, c.x);
-		else outp.arc = atan2(c.y, c.x) + 2.0*GCTL_Pi;
-		return outp;
-	}
-
-    template <typename T>
-	gctl::point2c<T> p2c(const gctl::point2p<T>& p)
-	{
-		point2c<T> outc;
-		outc.x = p.rad * cos(p.arc);
-		outc.y = p.rad * sin(p.arc);
-		return outc;
-	}
-
-    template <typename T>
-	gctl::point3s<T> c2s(const gctl::point3c<T>& c)
-	{
-		point3s<T> outs;
-		outs.rad = c.module();
-		if (outs.rad <= GCTL_ZERO) //点距离原点极近 将点置于原点
-		{
-			throw runtime_error("The point is at the origin. From point3c::c2s()");
-		}
-
-		outs.lat = 90.0 - acos(c.z/outs.rad)*180.0/GCTL_Pi;
-		outs.lon = atan2(c.y, c.x)*180.0/GCTL_Pi;
-		return outs;
-	}
-
-    template <typename T>
-	gctl::point3c<T> s2c(const gctl::point3s<T>& s)
-	{
-		//point3c<T> outc;
-		//outc.x = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*cos((2.0 + lon/180.0)*GCTL_Pi);
-		//outc.y = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*sin((2.0 + lon/180.0)*GCTL_Pi);
-		//outc.z = rad*cos((0.5 - lat/180.0)*GCTL_Pi);
-		//return outc;
-
-		point3c<T> v;
-		v.x = s.rad*cos(GCTL_Pi*s.lat/180.0)*cos(GCTL_Pi*s.lon/180.0);
-		v.y = s.rad*cos(GCTL_Pi*s.lat/180.0)*sin(GCTL_Pi*s.lon/180.0);
-		v.z = s.rad*sin(GCTL_Pi*s.lat/180.0);
-		return v;
-	}
-
     // this variable is only valid in this source file to control
     // the output of a 2D or 3D vertex.
     static char vertex_delimiter = ' ';
@@ -219,6 +168,6 @@ namespace gctl
         os >> a.id; a.in_loc(os);
         return os;
     }
-};
+}
 
 #endif // _GCTL_VERTEX_H

lib/graphic/cliplot.cpp

@@ -1,283 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#include "cliplot.h"
-
-gctl::cliplot::cliplot(size_t width, size_t height, double xmin, double xmax, double ymin, double ymax)
-{
-    new_screen_ = false;
-    digs_ = 4;
-    xl_num_ = 5;
-    yl_num_ = 5;
-    wname_ = "X";
-    hname_ = "Y";
-
-    struct winsize w;
-    ioctl(STDOUT_FILENO, TIOCGWINSZ, &w);
-    width += digs_ + 1;
-    width_ = w.ws_col<width?w.ws_col:width;
-    height_ =  w.ws_row<height?w.ws_row:height;
-    w0_ = digs_ + 1;
-    h0_ = height_ - 2;
-
-    xmin_ = xmin;
-    xmax_ = xmax;
-    ymin_ = ymin;
-    ymax_ = ymax;
-    dx_ = (xmax_ - xmin_)/(width_ - w0_);
-    dy_ = (ymax_ - ymin_)/height_;
-    sym_.resize(width_*height_);
-    att_.resize(width_*height_);
-    clear();
-    return;
-}
-
-gctl::cliplot::~cliplot(){}
-
-void gctl::cliplot::clear()
-{
-    std::fill(sym_.begin(), sym_.end(), ' ');
-    std::fill(att_.begin(), att_.end(), "");
-    return;
-}
-
-/**
- * @brief 设置画布字符属性
- * 
- * @param w 字符列位置
- * @param h 字符行位置
- * @param sym 字符
- * @param att 字符属性
- */
-void gctl::cliplot::set(size_t w, size_t h, char sym, const std::string &att)
-{
-    if (w >= width_ || h >= height_) return;
-    sym_[h*width_ + w] = sym;
-    att_[h*width_ + w] = att;
-    return;
-}
-
-void gctl::cliplot::set_axis(int xt_num, int yt_num)
-{
-    xl_num_ = xt_num;
-    yl_num_ = yt_num;
-    return;
-}
-
-void gctl::cliplot::set_digs(int digs)
-{
-    int old_digs = digs_;
-    digs_ = digs;
-
-    width_ += digs_ - old_digs;
-    w0_ += digs_ - old_digs;
-    sym_.resize(width_*height_);
-    att_.resize(width_*height_);
-    clear();
-    return;
-}
-
-void gctl::cliplot::set_new_screen(bool new_screen)
-{
-    new_screen_ = new_screen;
-    return;
-}
-
-void gctl::cliplot::set_wname(const std::string &wname)
-{
-    wname_ = wname;
-    return;
-}
-
-void gctl::cliplot::set_hname(const std::string &hname)
-{
-    hname_ = hname;
-    return;
-}
-
-void gctl::cliplot::plot_line(double x1, double x2, double y1, double y2, char s, const std::string &t)
-{
-    if (x1 > x2)
-    {
-        std::swap(x1, x2);
-        std::swap(y1, y2);
-    }
-
-    for (size_t w = w0_; w < width_; w++)
-    {
-        double x = xmin_ + (w - w0_)*(xmax_ - xmin_)/(width_ - 1 - w0_);
-        if (x + 0.5*dx_ >= x1 && x - 0.5*dx_ <= x2)
-        {
-            double y;
-            if (x2 - x1 < dx_) y = y1;
-            else y = y1 + (y2 - y1)*(x - 0.25*dx_ - x1)/(x2 - x1);
-            int h1 = round(h0_ - 1 - (y - ymin_)*(h0_ - 1)/(ymax_ - ymin_));
-
-            if (x2 - x1 < dx_) y = y2;
-            else y = y1 + (y2 - y1)*(x + 0.25*dx_ - x1)/(x2 - x1);
-            int h2 = round(h0_ - 1 - (y - ymin_)*(h0_ - 1)/(ymax_ - ymin_));
-
-            if (h1 > h2) std::swap(h1, h2);
-            for (size_t h = h1; h <= h2; h++)
-                set(w, h, s, t);
-        }
-    }
-    return;
-}
-
-void gctl::cliplot::plot_data(const std::vector<double> &x, const std::vector<double> &y, char s, const std::string &t)
-{
-    if (x.size() != y.size()) return;
-
-    for (size_t i = 0; i < x.size() - 1; i++)
-    {
-        plot_line(x[i], x[i+1], y[i], y[i+1], s, t);
-    }
-    return;
-}
-
-void gctl::cliplot::display(std::ostream &os)
-{
-    if (new_screen_)
-    {
-        os << GCTL_CLEARALL;
-        GCTL_MOVEUP(os, height_);
-    }
-
-    plot_axis();
-
-    for (size_t i = 0; i < height_; i++)
-    {
-        for (size_t j = 0; j < width_; j++)
-        {
-            os << att_[i*width_ + j] << sym_[i*width_ + j] << GCTL_RESET;
-        }
-        os << std::endl;
-    }
-    return;
-}
-
-std::string gctl::cliplot::axis_label(double num, int digs, int &odr)
-{
-    std::stringstream ss;
-    if (floor(log10(abs(num))) >= digs)
-    {
-        ss << std::scientific << std::setprecision(digs) << num;
-        odr = floor(log10(abs(num)));
-        if (num < 0) odr *= -1;
-    }
-    else
-    {
-        ss << std::fixed << std::setprecision(digs) << num;
-        odr = 0;
-    }
-    return ss.str();
-}
-
-void gctl::cliplot::plot_axis()
-{
-    for (size_t i = w0_; i < width_; i++)
-        set(i, h0_, '-', GCTL_BOLD);
-
-    std::vector<int> xpos(xl_num_);
-    int xtick = (width_ - w0_)/(xl_num_ - 1);
-    xpos[0] = w0_;
-    xpos[xl_num_ - 1] = width_ - 1;
-    for (size_t i = 1; i <= xl_num_ - 2; i++)
-        xpos[i] = w0_ + i*xtick;
-
-    int h_digs = digs_/2;
-    int odr;
-    for (size_t i = 0; i < xl_num_; i++)
-    {
-        set(xpos[i], h0_, '|', GCTL_BOLD);
-        std::string xlabel = axis_label(xmin_ + (xpos[i] - w0_)*(xmax_ - xmin_)/(width_ - 1 - w0_), digs_, odr);
-
-        for (int c = -h_digs; c <= h_digs; c++)
-        {
-            set(xpos[i] + c, h0_ + 1, xlabel[c + h_digs], GCTL_BOLD);
-        }
-
-        if (odr > 0)
-        {
-            set(xpos[i] + h_digs + 1, h0_ + 1, 'e', GCTL_BOLD);
-            set(xpos[i] + h_digs + 3, h0_ + 1, std::to_string(odr)[0], GCTL_BOLD);
-        }
-        else if (odr < 0)
-        {
-            set(xpos[i] + h_digs + 1, h0_ + 1, 'e', GCTL_BOLD);
-            set(xpos[i] + h_digs + 2, h0_ + 1, std::to_string(odr)[0], GCTL_BOLD);
-            set(xpos[i] + h_digs + 3, h0_ + 1, std::to_string(odr)[1], GCTL_BOLD);
-        }
-    }
-
-    for (size_t j = 0; j < h0_; j++)
-    {
-        set(w0_, j, '|', GCTL_BOLD);
-    }
-
-    std::vector<int> ypos(yl_num_);
-    int ytick = h0_/(yl_num_ - 1);
-    ypos[0] = 0;
-    ypos[yl_num_ - 1] = h0_ - 1;
-    for (size_t i = 1; i <= yl_num_ - 2; i++)
-        ypos[i] = i*ytick;
-
-    for (size_t i = 0; i < yl_num_; i++)
-    {
-        set(w0_, ypos[i], '-', GCTL_BOLD);
-        std::string ylabel = axis_label(ymax_ - ypos[i]*(ymax_ - ymin_)/(h0_ - 1), digs_, odr);
-
-        for (size_t c = 0; c < digs_; c++)
-        {
-            set(c, ypos[i], ylabel[c], GCTL_BOLD);
-        }
-
-        if (odr > 0)
-        {
-            set(digs_ + 3, ypos[i], 'e', GCTL_BOLD);
-            set(digs_ + 4, ypos[i], std::to_string(odr)[0], GCTL_BOLD);
-        }
-        else if (odr < 0)
-        {
-            set(digs_ + 3, ypos[i], 'e', GCTL_BOLD);
-            set(digs_ + 4, ypos[i], std::to_string(odr)[0], GCTL_BOLD);
-            set(digs_ + 5, ypos[i], std::to_string(odr)[1], GCTL_BOLD);
-        }
-    }
-
-    for (size_t i = 0; i < wname_.size(); i++)
-    {
-        set(width_ - wname_.size() + i, h0_ - 1, wname_[i], GCTL_BOLD);
-    }
-
-    for (size_t i = 0; i < hname_.size(); i++)
-    {
-        set(w0_ + 6 + i, 0, hname_[i], GCTL_BOLD);
-    }
-    return;
-}

lib/graphic/cliplot.h

@@ -1,206 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#ifndef _GCTL_CLIPLOT_H
-#define _GCTL_CLIPLOT_H
-
-#include <sys/ioctl.h>
-#include <unistd.h>
-#include <iostream>
-#include <sstream>
-#include <iomanip>
-#include <random>
-#include <algorithm>
-
-#include "../io/term_io.h"
-
-namespace gctl
-{
-    class cliplot
-    {
-    public:
-        /**
-         * @brief 构建一个新的画布对象
-         * 
-         * @param width 画布宽度（大于终端宽度则取终端宽度）
-         * @param height 画布高度（大于终端高度则取终端高度）
-         * @param xmin x轴最小值
-         * @param xmax x轴最大值
-         * @param ymin y轴最小值
-         * @param ymax y轴最大值
-         */
-        cliplot(size_t width, size_t height, double xmin, double xmax, double ymin, double ymax);
-
-        virtual ~cliplot(); ///< 析构函数
-
-        /**
-         * @brief 清除画布数据
-         */
-        void clear();
-
-        /**
-         * @brief 设置画布字符属性
-         * 
-         * @param w 字符列位置
-         * @param h 字符行位置
-         * @param sym 字符
-         * @param att 字符属性
-         */
-        void set(size_t w, size_t h, char sym, const std::string &att);
-
-        /**
-         * @brief 设置坐标轴刻度数量
-         * 
-         * @param xt_num x轴刻度数量
-         * @param yt_num y轴刻度数量
-         */
-        void set_axis(int xt_num, int yt_num);
-
-        /**
-         * @brief 设置坐标轴刻度数字位数
-         *
-         * @param digs 数字位数
-         */
-        void set_digs(int digs);
-
-        /**
-         * @brief 设置是否使用新的屏幕显示画布
-         *
-         * @param new_screen 是否使用新的屏幕显示画布
-         */
-        void set_new_screen(bool new_screen);
-
-        /**
-         * @brief 设置坐标轴名称
-         *
-         * @param wname x轴名称
-         */
-        void set_wname(const std::string &wname);
-
-        /**
-         * @brief 设置坐标轴名称
-         *
-         * @param hname y轴名称
-         */
-        void set_hname(const std::string &hname);
-
-        /**
-         * @brief 绘制直线
-         *
-         * @param x1 直线起点x坐标
-         * @param x2 直线终点x坐标
-         * @param y1 直线起点y坐标
-         * @param y2 直线终点y坐标
-         * @param s 直线字符
-         * @param t 直线字符属性
-         */
-        void plot_line(double x1, double x2, double y1, double y2, char s = '.', const std::string &t = GCTL_BOLDGREEN);
-
-        /**
-         * @brief 绘制函数
-         *
-         * @param func 函数
-         * @param s 函数字符
-         * @param t 函数字符属性
-         */
-        template <typename FuncOp>
-        void plot_func(FuncOp func, char s = '.', const std::string &t = GCTL_BOLDGREEN);
-
-        /**
-         * @brief 绘制数据
-         *
-         * @param x x坐标
-         * @param y y坐标
-         * @param s 数据点字符
-         * @param t 数据点字符属性
-         */
-        void plot_data(const std::vector<double> &x, const std::vector<double> &y, char s = '.', const std::string &t = GCTL_BOLDGREEN);
-
-        /**
-         * @brief 显示画布
-         *
-         * @param os 输出流
-         */
-        void display(std::ostream &os = std::cout);
-
-    private:
-        /**
-         * @brief 计算坐标轴刻度标签
-         *
-         * @param num 数值
-         * @param digs 数字位数
-         * @param odr 数值的10的幂次
-         * @return std::string 坐标轴刻度标签
-         */
-        std::string axis_label(double num, int digs, int &odr);
-
-        /**
-         * @brief 绘制坐标轴
-         */
-        void plot_axis();
-
-    private:
-        bool new_screen_;
-        int digs_; // 数字位数（包括小数点）
-        int xl_num_; // x轴刻度数量
-        int yl_num_; // y轴刻度数量
-        double xmin_; // x轴最小值
-        double xmax_; // x轴最大值
-        double ymin_; // y轴最小值
-        double ymax_; // y轴最大值
-        double dx_, dy_; // x轴与y轴的刻度间隔
-        size_t w0_, h0_; // 原点的行与列位置
-        size_t width_; // 画布宽度
-        size_t height_; // 画布高度
-        std::vector<char> sym_; // 画布字符
-        std::vector<std::string> att_; // 画布字符属性
-        std::string wname_;
-        std::string hname_;
-    };
-
-    template <typename FuncOp>
-    void gctl::cliplot::plot_func(FuncOp func, char s, const std::string &t)
-    {
-        for (size_t w = w0_; w < width_; w++)
-        {
-            double y, x = xmin_ + (w - w0_)*(xmax_ - xmin_)/(width_ - 1 - w0_);
-
-            y = func(x - 0.25*dx_);
-            int h1 = round(h0_ - 1 - (y - ymin_)*(h0_ - 1)/(ymax_ - ymin_));
-
-            y = func(x + 0.25*dx_);
-            int h2 = round(h0_ - 1 - (y - ymin_)*(h0_ - 1)/(ymax_ - ymin_));
-
-            if (h1 > h2) std::swap(h1, h2);
-            for (size_t h = h1; h <= h2; h++)
-                set(w, h, s, t);
-        }
-        return;
-    }
-};
-
-#endif // _GCTL_CLIPLOT_H

lib/graphic/gmt.cpp

@@ -1,148 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#include "gmt.h"
-
-gctl::gmt::gmt()
-{
-    seesion_ok_ = false;
-    GMT_API_ = nullptr;
-    session_name_ = "";
-    buffer_.clear();
-    tmpfiles_.clear();
-}
-
-gctl::gmt::~gmt()
-{
-    if (GMT_API_ != nullptr)
-    {
-        GMT_Destroy_Session(GMT_API_);
-        GMT_API_ = nullptr;
-        seesion_ok_ = false;
-        session_name_ = "";
-    }
-
-    if (!buffer_.empty()) buffer_.clear();
-    if (!tmpfiles_.empty()) tmpfiles_.clear();
-}
-
-void gctl::gmt::begin_session(const std::string& pic_name, 
-    const std::string& pic_type, const std::string& name)
-{
-    if (!buffer_.empty()) buffer_.clear();
-    session_name_ = pic_name;
-
-    GMT_API_ = GMT_Create_Session(name.c_str(), 2U, 0, nullptr);
-    if (GMT_API_ == nullptr) throw std::runtime_error("[gctl::gmt] Fail to begin GMT session.");
-
-    std::string cmd = pic_name + " " + pic_type;
-    if (GMT_Call_Module(GMT_API_, "begin", GMT_MODULE_CMD, (char*) cmd.c_str()) != GMT_NOERROR)
-        throw std::runtime_error("[gctl::gmt] Fail to execute GMT command.");
-    else
-    {
-        cmd = "begin " + cmd;
-        buffer_.push_back(cmd);
-    }
-
-    seesion_ok_ = true;
-    return;
-}
-
-void gctl::gmt::end_session(bool show, bool clear_tmp)
-{
-    std::string cmd = "";
-    if (show) cmd = "show";
-
-    if (GMT_Call_Module(GMT_API_, "end", GMT_MODULE_CMD, (char*) cmd.c_str()) != GMT_NOERROR)
-        throw std::runtime_error("[gctl::gmt] Fail to end GMT session.");
-    else
-    {
-        cmd = "end " + cmd;
-        buffer_.push_back(cmd);
-    }
-
-    if (clear_tmp)
-    {
-        for (auto& file : tmpfiles_) remove(file.c_str());
-    }
-
-    if (GMT_API_ != nullptr)
-    {
-        GMT_Destroy_Session(GMT_API_);
-        GMT_API_ = nullptr;
-    }
-
-    seesion_ok_ = false;
-    session_name_ = "";
-    return;
-}
-
-void gctl::gmt::save_session(const std::string& filename)
-{
-    time_t now = time(0);
-	char* dt = ctime(&now);
-
-    std::ofstream outfile;
-    open_outfile(outfile, filename, ".sh");
-
-    outfile << "#!/bin/bash" << std::endl;
-    outfile << "# This script is saved by gctl::gmt on " << dt;
-    outfile << "# For more information please connect yizhang-geo@zju.edu.cn" << std::endl;
-
-    for (auto& line : buffer_) outfile << "gmt " << line << std::endl;
-    outfile.close();
-
-    buffer_.clear();
-    return;
-}
-
-void gctl::gmt::call_module(const std::string& module, const std::string& cmd)
-{
-    if (GMT_Call_Module(GMT_API_, module.c_str(), GMT_MODULE_CMD, (char*) cmd.c_str()) != GMT_NOERROR)
-        throw std::runtime_error("[gctl::gmt] Fail to execute GMT command.");
-    else
-    {
-        std::string cmds = module + " " + cmd;
-        buffer_.push_back(cmds);
-    }
-}
-
-std::string gctl::gmt::create_grid(const array<double> &data, int xnum, int ynum, 
-    double xmin, double dx, double ymin, double dy)
-{
-    if (seesion_ok_ == false)
-        throw std::runtime_error("[gctl::gmt] GMT session is not started.");
-
-    if (xnum*ynum != data.size())
-        throw std::runtime_error("[gctl::gmt] The size of data does not match the grid size.");
-
-    std::string gridname = session_name_ + "_grid_" + std::to_string(tmpfiles_.size() + 1) + ".nc";
-    gctl::save_netcdf_grid(gridname, data, xnum, ynum, xmin, dx, ymin, dy);
-
-    tmpfiles_.push_back(gridname);
-    return gridname;
-}

lib/graphic/gmt.h

@@ -1,119 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#ifndef _GCTL_GMT_H
-#define _GCTL_GMT_H
-
-#include <vector>
-#include <exception>
-#include <stdexcept>
-#include <ctime>
-
-#include "gmt/gmt.h"
-#include "../io/file_io.h"
-#include "../io/netcdf_io.h"
-
-#ifndef GMT_VF_LEN
-#define GMT_VF_LEN 16
-#endif
-
-namespace gctl
-{
-    /**
-     * @brief 这是一个gmt命令管道库
-     */
-    class gmt
-    {
-    public:
-        /**
-         * @brief Construct a new gmt object
-         */
-        gmt();
-        // destructor
-        virtual ~gmt();
-
-        /**
-         * @brief Start a session.
-         * 
-         * @param pic_name Output picture name.
-         * @param pic_type Output picture type. Multiple types can be specified by comma.
-         * @param name Session name.
-         */
-        void begin_session(const std::string& pic_name, const std::string& pic_type = "eps", 
-            const std::string& name = "gctl_gmt_session");
-
-        /**
-         * @brief Destroy a session object
-         * 
-         * @param show Whether to show the result after finish the plotting.
-         * @param clear_tmp Whether to clear the temporary files.
-         */
-        void end_session(bool show = false, bool clear_tmp = true);
-
-        /**
-         * @brief Save the session to a file.
-         *
-         * @param filename File name.
-         */
-        void save_session(const std::string& filename);
-
-        /**
-         * @brief Call a GMT module.
-         * 
-         * @param module Name of the module.
-         * @param cmd Module command.
-         */
-        void call_module(const std::string& module, const std::string& cmd);
-
-        /**
-         * @brief Create a temporary grid file that could be used for plotting.
-         * 
-         * @param data Data array.
-         * @param xnum Number of data points in x direction.
-         * @param ynum Number of data points in y direction.
-         * @param xmin minimal x coordinate of the data.
-         * @param dx increment of x coordinate of the data.
-         * @param ymin minimal y coordinate of the data.
-         * @param dy increment of y coordinate of the data.
-         * @return Name of the virtual grid which could be used for plotting.
-         */
-        std::string create_grid(const array<double> &data, int xnum, int ynum, 
-            double xmin, double dx, double ymin, double dy);
-
-    private:
-        gmt(gmt const&) = delete;
-        void operator=(gmt const&) = delete;
-
-        bool seesion_ok_;
-        void* GMT_API_;
-        std::string session_name_;
-        std::vector<std::string> buffer_;
-        std::vector<std::string> tmpfiles_;
-    };
-};
-
-#endif // _GCTL_GMT_H

lib/graphic/gnuplot.cpp

@@ -1,135 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#include "gnuplot.h"
-
-gctl::gnuplot::gnuplot(bool persist)
-{
-    _2buffer = false;
-    pipe = popen(persist ? "gnuplot -persist" : "gnuplot", "w");
-    if (!pipe) throw std::runtime_error("gnuplot pipe failed to open. Try install gnuplot first.");
-}
-
-gctl::gnuplot::~gnuplot()
-{
-    if (pipe) pclose(pipe);
-}
-
-void gctl::gnuplot::to_buffer()
-{
-    _2buffer = true;
-    return;
-}
-
-void gctl::gnuplot::send(const std::string& cmd, bool next_block)
-{
-    if (!pipe) return;
-
-    if (!_2buffer) fputs((cmd + "\n").c_str(), pipe);
-    else
-    {
-        if (next_block) buffer.push_back(cmd + "\n\n");
-        else buffer.push_back(cmd + "\n");
-    }
-    return;
-}
-
-void gctl::gnuplot::send_buffer(size_t repeat)
-{
-    if (!pipe) return;
-    for (size_t i = 0; i < repeat; i++)
-    {
-        for (auto& line : buffer) fputs(line.c_str(), pipe);
-        fputs("\n", pipe);
-    }
-
-    fflush(pipe);
-    buffer.clear();
-    return;
-}
-
-void gctl::gnuplot::send_data(std::string name, const std::vector<std::vector<double> > &data, bool next_block)
-{
-    if (!pipe) return;
-
-    if (!_2buffer)
-    {
-        fputs((name + " << EOD\n").c_str(), pipe);
-
-        int cnum = data.size();
-        int dnum = data[0].size();
-        std::string line;
-        for (size_t i = 0; i < dnum; i++)
-        {
-            line = std::to_string(data[0][i]);
-            for (size_t j = 1; j < cnum; j++)
-            {
-                line += " " + std::to_string(data[j][i]); 
-            }
-            fputs((line + "\n").c_str(), pipe);
-        }
-
-        fputs("EOD\n", pipe);
-    }
-    else
-    {
-        buffer.push_back(name + " << EOD\n");
-
-        int cnum = data.size();
-        int dnum = data[0].size();
-        std::string line;
-        for (size_t i = 0; i < dnum; i++)
-        {
-            line = std::to_string(data[0][i]);
-            for (size_t j = 1; j < cnum; j++)
-            {
-                line += " " + std::to_string(data[j][i]); 
-            }
-            buffer.push_back(line + "\n");
-        }
-
-        buffer.push_back("EOD\n");
-        if (next_block) buffer.push_back("\n");
-    }
-    return;
-}
-
-void gctl::gnuplot::save_buffer(const std::string& filename)
-{
-    time_t now = time(0);
-	char* dt = ctime(&now);
-
-    std::ofstream outfile;
-    open_outfile(outfile, filename, ".gp");
-
-    outfile << "# This script is saved by gctl::gnuplot on " << dt;
-    outfile << "# For more information please connect yizhang-geo@zju.edu.cn" << std::endl;
-
-    for (auto& line : buffer) outfile << line;
-    outfile.close();
-    return;
-}

lib/graphic/gnuplot.h

@@ -1,103 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#ifndef _GCTL_GNUPLOT_H
-#define _GCTL_GNUPLOT_H
-
-#include <vector>
-#include <exception>
-#include <stdexcept>
-#include <ctime>
-
-#include "../io/file_io.h"
-
-namespace gctl
-{
-    /**
-     * @brief 这是一个gnuplot命令管道库
-     */
-    class gnuplot
-    {
-    public:
-        /**
-         * @brief Construct a new gnuplot object
-         * 
-         * @param persist Gnuplot will hold the window after plot if true.
-         */
-        gnuplot(bool persist = true);
-
-        // destructor
-        virtual ~gnuplot();
-
-        /**
-         * @brief Set this to send all lines to buffer.
-         */
-        void to_buffer();
-
-        /**
-         * @brief Send a line to gnuplot and execute it. If to_buffer is called before, then the line will be stored in the buffer.
-         * 
-         * @param cmd command line.
-         * @param next_block If true, an empty line will be added to the buffer after the cmd.
-         */
-        void send(const std::string& cmd, bool next_block = false);
-
-        /**
-         * @brief Send a data block
-         * 
-         * @param name Name of the data block which could be used later for ploting.
-         * @param data data block. Note that each vector object represents a data column.
-         * @param next_block If true, an empty line will be added to the buffer after the cmd.
-         */
-        void send_data(std::string name, const std::vector<std::vector<double> > &data, bool next_block = false);
-
-        /**
-         * @brief Send the buffer to gnuplot and execute it at once. Then clear the buffer.
-         * 
-         * @param repeat How many times to repeat the buffer.
-         */
-        void send_buffer(size_t repeat = 1);
-
-        /**
-         * @brief Save the buffer to a file. This will not clear the buffer.
-         *
-         * @param filename File name.
-         */
-        void save_buffer(const std::string& filename);
-
-    private:
-        gnuplot(gnuplot const&) = delete;
-        void operator=(gnuplot const&) = delete;
-
-        FILE* pipe;
-        std::vector<std::string> buffer;
-        std::vector<std::string> data_lines;
-        bool _2buffer;
-    };
-};
-
-#endif // _GCTL_GNUPLOT_H

lib/io.h

@@ -28,21 +28,15 @@
 #ifndef _GCTL_IO_H
 #define _GCTL_IO_H
 
-#include "io/term_io.h"
+#include "io/mesh_io.h"
-#include "io/file_io.h"
 #include "io/native_io.h"
-
-#include "io/surfer_io.h"
-#include "io/netcdf_io.h"
-
-#include "io/off_io.h"
-#include "io/stl_io.h"
-#include "io/ply_io.h"
-#include "io/tetgen_io.h"
-#include "io/triangle_io.h"
-#include "io/gmsh_io.h"
-
 #include "io/text_io.h"
 #include "io/dsv_io.h"
+#include "io/gmsh_io.h"
+#include "io/surfer_io.h"
+#include "io/tetgen_io.h"
+#include "io/triangle_io.h"
+#include "io/off_io.h"
+#include "io/netcdf_io.h"
 
 #endif // _GCTL_IO_H

lib/io/dsv_io.cpp

@@ -29,6 +29,7 @@
 
 gctl::dsv_io::dsv_io()
 {
+    file_ = "";
     att_sym_ = '#';
 	tag_sym_ = '!';
 	deli_sym_ = ' ';
@@ -42,8 +43,9 @@ gctl::dsv_io::~dsv_io()
     clear();
 }
 
-gctl::dsv_io::dsv_io(std::string filename, std::string file_exten, int t)
+gctl::dsv_io::dsv_io(std::string filename, std::string file_exten, table_headtype_e t)
 {
+    file_ = "";
     att_sym_ = '#';
 	tag_sym_ = '!';
 	deli_sym_ = ' ';
@@ -57,6 +59,7 @@ gctl::dsv_io::dsv_io(std::string filename, std::string file_exten, int t)
 
 void gctl::dsv_io::clear()
 {
+    file_ = "";
     att_sym_ = '#';
 	tag_sym_ = '!';
 	deli_sym_ = ' ';
@@ -71,17 +74,27 @@ void gctl::dsv_io::clear()
 	return;
 }
 
-std::vector<std::string> gctl::dsv_io::row_names()
+void gctl::dsv_io::get_row_names(std::vector<std::string> &names)
 {
-    _1s_vector names(row_num_);
+    names.resize(row_num_);
     for (size_t i = 1; i <= row_num_; i++)
     {
         names[i - 1] = table_[i][0].str_;
     }
-    return names;
+    return;
 }
 
-void gctl::dsv_io::row_names(const std::vector<std::string> &names, const std::vector<int> &idx, std::string corner_name)
+void gctl::dsv_io::get_column_names(std::vector<std::string> &names)
+{
+    names.resize(col_num_);
+    for (size_t i = 1; i <= col_num_; i++)
+    {
+        names[i - 1] = table_[0][i].str_;
+    }
+    return;
+}
+
+void gctl::dsv_io::set_row_names(const std::vector<std::string> &names, const std::vector<int> &idx, std::string corner_name)
 {
     if (!idx.empty())
     {
@@ -104,17 +117,7 @@ void gctl::dsv_io::row_names(const std::vector<std::string> &names, const std::v
     return;
 }
 
-std::vector<std::string> gctl::dsv_io::column_names()
+void gctl::dsv_io::set_column_names(const std::vector<std::string> &names, const std::vector<int> &idx)
-{
-    _1s_vector names(col_num_);
-    for (size_t i = 1; i <= col_num_; i++)
-    {
-        names[i - 1] = table_[0][i].str_;
-    }
-    return names;
-}
-
-void gctl::dsv_io::column_names(const std::vector<std::string> &names, const std::vector<int> &idx)
 {
     if (!idx.empty())
     {
@@ -135,11 +138,11 @@ void gctl::dsv_io::column_names(const std::vector<std::string> &names, const std
     return;
 }
 
-void gctl::dsv_io::row_type(table_cell_type t, int idx)
+void gctl::dsv_io::set_row_type(cell_type_e t, int idx)
 {
     if (idx > row_num_ || idx <= 0)
     {
-        throw std::runtime_error("[gctl::dsv_io] Invalid row index.");
+        throw std::runtime_error("[gctl::dsv_io] Invalid column index.");
     }
 
     for (size_t i = 0; i <= col_num_; i++)
@@ -148,27 +151,12 @@ void gctl::dsv_io::row_type(table_cell_type t, int idx)
     }
 }
 
-void gctl::dsv_io::row_type(table_cell_type t, std::string name)
+void gctl::dsv_io::set_row_type(cell_type_e t, std::string name)
 {
-    row_type(t, name_index(name, true));
+    set_row_type(t, name_index(name, true));
 }
 
-gctl::table_cell_type gctl::dsv_io::row_type(int idx)
+void gctl::dsv_io::set_column_type(cell_type_e t, int idx)
-{
-    if (idx > row_num_ || idx <= 0)
-    {
-        throw std::runtime_error("[gctl::dsv_io] Invalid row index.");
-    }
-
-    return table_[idx][0].type_;
-}
-
-gctl::table_cell_type gctl::dsv_io::row_type(std::string name)
-{
-    return row_type(name_index(name, true));
-}
-
-void gctl::dsv_io::column_type(table_cell_type t, int idx)
 {
     if (idx > col_num_ || idx <= 0)
     {
@@ -182,53 +170,39 @@ void gctl::dsv_io::column_type(table_cell_type t, int idx)
     return;
 }
 
-void gctl::dsv_io::column_type(table_cell_type t, std::string name)
+void gctl::dsv_io::set_column_type(cell_type_e t, std::string name)
 {
-    column_type(t, name_index(name, false));
+    set_column_type(t, name_index(name, false));
 }
 
-gctl::table_cell_type gctl::dsv_io::column_type(int idx)
+void gctl::dsv_io::load_text(std::string filename, std::string file_exten, table_headtype_e t)
-{
-    if (idx > col_num_ || idx <= 0)
-    {
-        throw std::runtime_error("[gctl::dsv_io] Invalid column index.");
-    }
-
-    return table_[0][idx].type_;
-}
-
-gctl::table_cell_type gctl::dsv_io::column_type(std::string name)
-{
-    return column_type(name_index(name, false));
-}
-
-void gctl::dsv_io::load_text(std::string filename, std::string file_exten, int t)
 {
     std::ifstream infile;
     open_infile(infile, filename, file_exten);
 
+    file_ = filename + file_exten;
+
     int h = 0;
     std::string tmp_line;
     std::vector<std::string> lines;
     while (std::getline(infile, tmp_line))
     {
         if (tmp_line.empty()) continue; // 跳过空行 空行不会并记入头信息计数中
-
+        else if (tmp_line[0] == att_sym_) // 注释行或者标记行 # #!
-        // 去掉空格
-        tmp_line.erase(0, tmp_line.find_first_not_of(" \t"));
-        tmp_line.erase(tmp_line.find_last_not_of(" \t") + 1);
-
-        if (tmp_line[0] == att_sym_) // 注释行或者标记行 # #!
         {
             if (tmp_line[1] == tag_sym_) // #!
             {
                 tmp_line = tmp_line.substr(2); // 去掉前两个字符
+                tmp_line.erase(0, tmp_line.find_first_not_of(" \t"));
+                tmp_line.erase(tmp_line.find_last_not_of(" \t") + 1);
                 tags_.push_back(tmp_line);
                 continue;
             }
-
+            
             // #
             tmp_line = tmp_line.substr(1); // 去掉第一个字符
+            tmp_line.erase(0, tmp_line.find_first_not_of(" \t"));
+            tmp_line.erase(tmp_line.find_last_not_of(" \t") + 1);
             annotates_.push_back(tmp_line);
         }
         else if (h < head_num_) //读入到头信息中
@@ -251,7 +225,7 @@ void gctl::dsv_io::load_text(std::string filename, std::string file_exten, int t
         
         // 动态调整列数
         cn = tmp_cols.size();
-        cn_max = GCTL_MAX(cn, cn_max);
+        cn_max = std::max(cn, cn_max);
 
         table_[i].resize(tmp_cols.size());
         for (size_t j = 0; j < tmp_cols.size(); j++)
@@ -286,7 +260,7 @@ void gctl::dsv_io::load_text(std::string filename, std::string file_exten, int t
         }
     }
 
-    if (t == ColHead) // 有列头 需要补齐空白的行头
+    if (t == ColumnHead) // 有列头 需要补齐空白的行头
     {
         row_num_ = table_.size() - 1;
         col_num_ = table_[0].size();
@@ -306,7 +280,7 @@ void gctl::dsv_io::load_text(std::string filename, std::string file_exten, int t
         table_[0].resize(col_num_ + 1);
     }
     
-    if ((t & RowHead) && (t & ColHead)) // 有行头和列头
+    if (t == BothHead) // 有行头和列头
     {
         row_num_ = table_.size() - 1;
         col_num_ = table_[0].size() - 1;
@@ -316,9 +290,9 @@ void gctl::dsv_io::load_text(std::string filename, std::string file_exten, int t
     return;
 }
 
-void gctl::dsv_io::load_csv(std::string filename, int t)
+void gctl::dsv_io::load_csv(std::string filename, table_headtype_e t)
 {
-    delimeter(',');
+    set_delimeter(',');
     load_text(filename, ".csv", t);
     return;
 }
@@ -343,30 +317,80 @@ void gctl::dsv_io::save_text(std::string filename, std::string file_exten)
         outfile << "# " << annotates_[i] << std::endl;
     }
 
+    // 探测是否有行头
+    bool col_st = 1;
     for (int i = 0; i <= row_num_; i++)
     {
-        for (size_t j = 0; j <= col_num_; j++)
+        if (table_[i][0].out_ok_ && table_[i][0].str_ != "")
         {
-            if (table_[i][j].out_ok_ && table_[i][j].str_!= "")
+            col_st = 0;
-            {
+            break;
-                outfile << table_[i][j].str_;
-                for (size_t k = j + 1; k <= col_num_; k++)
-                {
-                    if (table_[i][k].out_ok_) outfile << deli_sym_ << table_[i][k].str_;
-                }
-                outfile << std::endl;
-                break;
-            }
         }
     }
 
+    for (int i = 0; i <= row_num_; i++)
+    {
+        // 单独处理第一列 即行头
+        outfile << table_[i][col_st].str_;
+        for (int j = col_st + 1; j <= col_num_; j++)
+        {
+            if (table_[i][j].out_ok_) outfile << deli_sym_ << table_[i][j].str_;
+        }
+        outfile << std::endl;
+    }
+    
+/*
+    // 单独处理第一行 即列头
+    bool line_st = false;
+    if (table_[0][0].out_ok_ && table_[0][0].str_ != "")
+    {
+        outfile << table_[0][0].str_;
+        line_st = true;
+    }
+
+    for (int j = 1; j <= col_num_; j++)
+    {
+        if (line_st && table_[0][j].out_ok_ && table_[0][j].str_ != "") outfile << deli_sym_ << table_[0][j].str_; // line started
+        else if (table_[0][j].out_ok_ && table_[0][j].str_ != "") // line not started
+        {
+            outfile << table_[0][j].str_;
+            line_st = true; // start line
+        }
+    }
+
+    if (line_st) outfile << std::endl;
+
+    // 处理余下的行
+    for (int i = 1; i <= row_num_; i++)
+    {
+        line_st = false;
+
+        // 单独处理第一列 即行头
+        if (table_[i][0].out_ok_ && table_[i][0].str_ != "")
+        {
+            outfile << table_[i][0].str_;
+            line_st = true;
+        }
+
+        for (int j = 1; j <= col_num_; j++)
+        {
+            if (line_st && table_[i][j].out_ok_) outfile << deli_sym_ << table_[i][j].str_; // line started
+            else if (table_[i][j].out_ok_) // line not started
+            {
+                outfile << table_[i][j].str_;
+                line_st = true; // start line
+            }
+        }
+        outfile << std::endl;
+    }
+*/
     outfile.close();
     return;
 }
 
 void gctl::dsv_io::save_csv(std::string filename)
 {
-    delimeter(',');
+    set_delimeter(',');
     save_text(filename, ".csv");
     return;
 }
@@ -387,167 +411,72 @@ void gctl::dsv_io::init_table(int row, int col)
     return;
 }
 
-gctl::dsv_io gctl::dsv_io::export_table(bool ignore_disabled)
+void gctl::dsv_io::info(table_headtype_e t)
 {
-    std::vector<std::string> str_line, row_names, col_names;
+    std::clog << "File: " << file_ << "\n------------\n";
-    std::vector<std::vector<std::string> > str_table;
+    std::clog << "Head(s): " << head_num_ << "\n";
+    std::clog << "Annotation(s): " << annotates_.size() << "\n";
+    std::clog << "Tag(s): " << tags_.size() << "\n";
+    std::clog << "------------\n";
 
-    std::string cor_name = table_[0][0].str_;
+    if (t == ColumnHead || t == BothHead)
-    for (size_t j = 1; j <= col_num_; j++)
     {
-        if (table_[0][j].out_ok_ || !ignore_disabled)
+        std::clog << "Columns:\n";
-            col_names.push_back(table_[0][j].str_);
-    }
-
-    for (size_t i = 1; i <= row_num_; i++)
-    {
-        if (table_[i][0].out_ok_ || !ignore_disabled)
-        {
-            str_line.clear();
-
-            for (size_t j = 1; j <= col_num_; j++)
-            {
-                if (table_[i][j].out_ok_ || !ignore_disabled)
-                    str_line.push_back(table_[i][j].str_);
-            }
-
-            str_table.push_back(str_line);
-            row_names.push_back(table_[i][0].str_);
-        }
-    }
-
-    dsv_io out_table;
-    out_table.init_table(str_table);
-    out_table.row_names(row_names, {}, cor_name);
-    out_table.column_names(col_names);
-
-    destroy_vector(row_names);
-    destroy_vector(col_names);
-    destroy_vector(str_line);
-    destroy_vector(str_table);
-    return out_table;
-}
-
-void gctl::dsv_io::info(int t, std::ostream &os)
-{
-    if (t & HeadInfo)
-    {
-        os << "Head(s): " << head_num_ << "\n";
-        for (size_t i = 0; i < heads_.size(); i++)
-        {
-            os << heads_[i] << "\n";
-        }
-        os << "------------\n";
-    }
-
-    if (t & AttInfo)
-    {
-        os << "Annotation(s): " << annotates_.size() << "\n";
-        for (size_t i = 0; i < annotates_.size(); i++)
-        {
-            os << annotates_[i] << "\n";
-        }
-        os << "------------\n";
-    }
-
-    if (t & TagInfo)
-    {
-        os << "Tag(s): " << tags_.size() << "\n";
-        for (size_t i = 0; i < tags_.size(); i++)
-        {
-            os << tags_[i] << "\n";
-        }
-        os << "------------\n";
-    }
-
-    if (t & ColInfo)
-    {
-        os << "Columns:\n";
         for (size_t i = 1; i <= col_num_; i++)
         {
             if (table_[0][i].str_ != "")
             {
-                os << table_[0][i].str_ << " | ";
+                std::clog << table_[0][i].str_ << " | ";
-                if (table_[0][i].out_ok_) os << "Enabled | ";
+                if (table_[1][i].type_ == String) std::clog << "String | ";
-                else os << "Disabled | ";
+                if (table_[1][i].type_ == Int) std::clog << "Int | ";
-                if (table_[1][i].type_ == String) os << "String | ";
+                if (table_[1][i].type_ == Float) std::clog << "Float | ";
-                if (table_[1][i].type_ == Int) os << "Int | ";
+                std::clog << table_[1][i].str_ << " -> " << table_[row_num_][i].str_;
-                if (table_[1][i].type_ == Float) os << "Float | ";
-                os << table_[1][i].str_ << " -> " << table_[row_num_][i].str_;
             }
             else
             {
-                os << "C" + std::to_string(i) << " | ";
+                std::clog << "C" + std::to_string(i) << " | ";
-                if (table_[0][i].out_ok_) os << "Enabled | ";
+                if (table_[1][i].type_ == String) std::clog << "String | ";
-                else os << "Disabled | ";
+                if (table_[1][i].type_ == Int) std::clog << "Int | ";
-                if (table_[1][i].type_ == String) os << "String | ";
+                if (table_[1][i].type_ == Float) std::clog << "Float | ";
-                if (table_[1][i].type_ == Int) os << "Int | ";
+                std::clog << table_[1][i].str_ << " -> " << table_[row_num_][i].str_;
-                if (table_[1][i].type_ == Float) os << "Float | ";
-                os << table_[1][i].str_ << " -> " << table_[row_num_][i].str_;
             }
-            os << std::endl;
+
+            if (!table_[0][i].out_ok_) std::clog << " (No output)";
+            std::clog << std::endl;
         }
-        os << "------------\n";
+        std::clog << "============\n";
     }
     
-    if (t & RowInfo)
+    if (t == RowHead || t == BothHead)
     {
-        os << "Rows:\n";
+        std::clog << "Rows:\n";
         for (size_t i = 1; i <= row_num_; i++)
         {
             if (table_[i][0].str_ != "")
             {
-                os << table_[i][0].str_ << " | ";
+                std::clog << table_[i][0].str_ << " | ";
-                if (table_[i][0].out_ok_) os << "Enabled | ";
+                if (table_[i][1].type_ == String) std::clog << "String | ";
-                else os << "Disabled | ";
+                if (table_[i][1].type_ == Int) std::clog << "Int | ";
-                if (table_[i][1].type_ == String) os << "String | ";
+                if (table_[i][1].type_ == Float) std::clog << "Float | ";
-                if (table_[i][1].type_ == Int) os << "Int | ";
+                std::clog << table_[i][1].str_ << " -> " << table_[i][col_num_].str_;
-                if (table_[i][1].type_ == Float) os << "Float | ";
-                os << table_[i][1].str_ << " -> " << table_[i][col_num_].str_;
             }
             else
             {
-                os << "R" + std::to_string(i) << " | ";
+                std::clog << "R" + std::to_string(i) << " | ";
-                if (table_[i][0].out_ok_) os << "Enabled | ";
+                if (table_[i][1].type_ == String) std::clog << "String | ";
-                else os << "Disabled | ";
+                if (table_[i][1].type_ == Int) std::clog << "Int | ";
-                if (table_[i][1].type_ == String) os << "String | ";
+                if (table_[i][1].type_ == Float) std::clog << "Float | ";
-                if (table_[i][1].type_ == Int) os << "Int | ";
+                std::clog << table_[i][1].str_ << " -> " << table_[i][col_num_].str_;
-                if (table_[i][1].type_ == Float) os << "Float | ";
-                os << table_[i][1].str_ << " -> " << table_[i][col_num_].str_;
             }
-            os << std::endl;
+
+            if (!table_[i][0].out_ok_) std::clog << " (No output)";
+            std::clog << std::endl;
         }
-        os << "------------\n";
+        std::clog << "============\n";
     }
     return;
 }
 
-void gctl::dsv_io::display()
-{
-    std::string line;
-    cli_viewer viewer;
-    for (int i = 0; i <= row_num_; i++)
-    {
-        for (size_t j = 0; j <= col_num_; j++)
-        {
-            if (table_[i][j].out_ok_ && table_[i][j].str_!= "")
-            {
-                line = table_[i][j].str_;
-                for (size_t k = j + 1; k <= col_num_; k++)
-                {
-                    if (table_[i][k].out_ok_) line += deli_sym_ + table_[i][k].str_;
-                }
-
-                viewer.addData(line);
-                break;
-            }
-        }
-    }
-
-    viewer.display();
-    return;
-}
-
 int gctl::dsv_io::name_index(std::string name, bool iter_row)
 {
     // 拾取行号或列号 格式为R<id>和C<id>
@@ -588,19 +517,6 @@ int gctl::dsv_io::name_index(std::string name, bool iter_row)
     }
 }
 
-void gctl::dsv_io::table_output(switch_type_e s)
-{
-    for (size_t i = 0; i <= row_num_; i++)
-    {
-        for (size_t j = 0; j <= col_num_; j++)  
-        {
-            if (s == Enable) table_[i][j].out_ok_ = true;
-            else table_[i][j].out_ok_ = false;
-        }
-    }
-    return;
-}
-
 void gctl::dsv_io::column_output(int idx, switch_type_e s)
 {
     if (idx > col_num_ || idx <= 0)
@@ -707,63 +623,162 @@ int gctl::dsv_io::add_row(std::string name, std::string id_name)
     return add_row(name, name_index(id_name, true));
 }
 
-void gctl::dsv_io::filter(std::string cnd_str, std::string cnd_tar, table_headtype_e thead)
+void gctl::dsv_io::filt_column(std::string cnd_str, std::string cnd_col, 
+	const std::vector<std::string> &out_col, dsv_io &out_table)
 {
-    int idx;
+    int idx = name_index(cnd_col);
-    if (thead == RowHead) idx = name_index(cnd_tar, true);
+    if (idx < 0) throw std::runtime_error("[gctl::dsv_io::] Invalid column index or name.");
-    else if (thead == ColHead) idx = name_index(cnd_tar);
-    else throw std::runtime_error("[gctl::dsv_io::filter] Invalid table head type.");
 
-    if (idx < 0) throw std::runtime_error("[gctl::dsv_io::filter] Invalid row/column index or name.");
+    array<int> odx;
+    bool out_row = false;
+    if (out_col.empty()) out_row = true;
+    else
+    {
+        odx.resize(out_col.size());
+        for (size_t i = 0; i < out_col.size(); i++)
+        {
+            odx[i] = name_index(out_col[i]);
+            if (odx[i] < 0) throw std::runtime_error("[gctl::dsv_io::] Invalid column index or name.");
+        }
+    }
 
     std::smatch ret;
     std::regex pat(cnd_str);
-    if (thead == RowHead) // cnd_tar是行头 此时为按列过滤
+    std::vector<std::string> str_line, row_names;
+    std::vector<std::vector<std::string> > str_table; 
+
+    for (size_t i = 1; i <= row_num_; i++)
     {
-        for (size_t i = 1; i <= col_num_; i++)
+        if (regex_search(table_[i][idx].str_, ret, pat))
         {
-            if (!regex_search(table_[idx][i].str_, ret, pat))
+            if (out_row)
             {
-                column_output(i, Disable);
+                str_line.clear();
+
+                for (size_t j = 1; j <= col_num_; j++)
+                {
+                    str_line.push_back(table_[i][j].str_);
+                }
+                
+                str_table.push_back(str_line);
             }
+            else
+            {
+                str_line.clear();
+                str_line.push_back(table_[i][idx].str_);
+                
+                for (size_t j = 0; j < odx.size(); j++)
+                {
+                    str_line.push_back(table_[i][odx[j]].str_);
+                }
+                
+                str_table.push_back(str_line);
+            }
+
+            row_names.push_back(table_[i][0].str_);
         }
     }
-    else // cnd_tar是列头 此时为按行过滤
+
+    out_table.init_table(str_table);
+
+    std::vector<std::string> io_col;
+    if (out_row)
     {
-        for (size_t i = 1; i <= row_num_; i++)
+        get_column_names(io_col);
+        out_table.cell(table_[0][0].str_, 0, 0);
+    }
+    else
+    {
+        io_col.push_back(cnd_col);
+
+        for (size_t j = 0; j < odx.size(); j++)
         {
-            if (!regex_search(table_[i][idx].str_, ret, pat))
+            io_col.push_back(out_col[j]);
-            {
-                row_output(i, Disable);
-            }
         }
     }
+
+    out_table.set_column_names(io_col);
+    out_table.set_row_names(row_names, {}, table_[0][0].str_);
+
+    destroy_vector(row_names);
+    destroy_vector(io_col);
+    destroy_vector(str_line);
+    destroy_vector(str_table);
     return;
 }
 
-void gctl::dsv_io::filter(linebool_func_t func, table_headtype_e thead)
+void gctl::dsv_io::filt_column(rowbool_func_t func, const std::vector<std::string> &out_col, dsv_io &out_table)
 {
-    if (thead == RowHead)
+    array<int> odx;
+    bool out_row = false;
+    if (out_col.empty()) out_row = true;
+    else
     {
-        for (size_t i = 1; i <= row_num_; i++)
+        odx.resize(out_col.size());
+        for (size_t i = 0; i < out_col.size(); i++)
         {
-            if (!func(table_[i])) row_output(i, Disable);
+            odx[i] = name_index(out_col[i]);
+            if (odx[i] < 0) throw std::runtime_error("[gctl::dsv_io::] Invalid column index or name.");
         }
     }
-    else if (thead == ColHead)
+
+    std::vector<std::string> str_line, row_names;
+    std::vector<std::vector<std::string> > str_table; 
+
+    for (size_t i = 1; i <= row_num_; i++)
     {
-        std::vector<table_cell> col_cell(row_num_);
+        if (func(table_[i]))
-        for (size_t i = 1; i <= col_num_; i++)
         {
-            for (size_t j = 1; j < row_num_; j++)
+            if (out_row)
             {
-                col_cell[j] = table_[j][i];
+                str_line.clear();
+
+                for (size_t j = 1; j <= col_num_; j++)
+                {
+                    str_line.push_back(table_[i][j].str_);
+                }
+                
+                str_table.push_back(str_line);
             }
-            
+            else
-            if (!func(col_cell)) column_output(i, Disable);
+            {
+                str_line.clear();
+                
+                for (size_t j = 0; j < odx.size(); j++)
+                {
+                    str_line.push_back(table_[i][odx[j]].str_);
+                }
+                
+                str_table.push_back(str_line);
+            }
+
+            row_names.push_back(table_[i][0].str_);
         }
     }
-    else throw std::runtime_error("[gctl::dsv_io::filter] Invalid table head type.");
+
+    out_table.init_table(str_table);
+
+    std::vector<std::string> io_col;
+    if (out_row)
+    {
+        get_column_names(io_col);
+        out_table.cell(table_[0][0].str_, 0, 0);
+    }
+    else
+    {
+        for (size_t j = 0; j < odx.size(); j++)
+        {
+            io_col.push_back(out_col[j]);
+        }
+    }
+
+    out_table.set_column_names(io_col);
+    out_table.set_row_names(row_names, {}, table_[0][0].str_);
+
+    destroy_vector(row_names);
+    destroy_vector(io_col);
+    destroy_vector(str_line);
+    destroy_vector(str_table);
     return;
 }
 
@@ -783,7 +798,7 @@ void gctl::dsv_io::cal_column(std::string expr_str, const std::vector<std::strin
             throw std::runtime_error("[gctl::dsv_io] Invalid column type for numerical calculating.");
         }
     }
-
+    
     exprtk::symbol_table<double> symbol_table;
 	array<double> var(col_list.size());
 
@@ -815,45 +830,41 @@ void gctl::dsv_io::cal_column(std::string expr_str, const std::vector<std::strin
     return;
 }
 
-void gctl::dsv_io::filter(std::string cnd_str, const std::vector<std::string> &cnd_tars, table_headtype_e thead)
+void gctl::dsv_io::filt_column(std::string cnd_str, const std::vector<std::string> &cnd_col, 
+	const std::vector<std::string> &out_col, dsv_io& out_table)
 {
-    array<int> idx(cnd_tars.size());
+    array<int> idx(cnd_col.size());
-    if (thead == RowHead)
+    for (size_t i = 0; i < cnd_col.size(); i++)
     {
-        for (size_t i = 0; i < cnd_tars.size(); i++)
+        idx[i] = name_index(cnd_col[i]);
+
+        if (idx[i] < 0) throw std::runtime_error("[gctl::dsv_io::] Invalid column index or name.");
+
+        if (table_[0][idx[i]].type_ != Int && table_[0][idx[i]].type_ != Float)
         {
-            idx[i] = name_index(cnd_tars[i], true);
+            throw std::runtime_error("[gctl::dsv_io] Invalid column type for numerical calculating.");
-
-            if (idx[i] <= 0 || idx[i] > row_num_) throw std::runtime_error("[gctl::dsv_io::filter] Invalid row index or name.");
-
-            if (table_[idx[i]][0].type_ != Int && table_[idx[i]][0].type_ != Float)
-            {
-                throw std::runtime_error("[gctl::dsv_io::filter] Invalid row type for numerical calculating.");
-            }
         }
     }
-    else if (thead == ColHead)
+
+    array<int> odx;
+    bool out_row = false;
+    if (out_col.empty()) out_row = true;
+    else
     {
-        for (size_t i = 0; i < cnd_tars.size(); i++)
+        odx.resize(out_col.size());
+        for (size_t i = 0; i < out_col.size(); i++)
         {
-            idx[i] = name_index(cnd_tars[i]);
+            odx[i] = name_index(out_col[i]);
-
+            if (odx[i] < 0) throw std::runtime_error("[gctl::dsv_io::] Invalid column index or name.");
-            if (idx[i] <= 0 || idx[i] > col_num_) throw std::runtime_error("[gctl::dsv_io::filter] Invalid column index or name.");
-
-            if (table_[0][idx[i]].type_ != Int && table_[0][idx[i]].type_ != Float)
-            {
-                throw std::runtime_error("[gctl::dsv_io::filter] Invalid column type for numerical calculating.");
-            }
         }
     }
-    else throw std::runtime_error("[gctl::dsv_io::filter] Invalid table head type.");
+    
-
     exprtk::symbol_table<double> symbol_table;
-	array<double> var(cnd_tars.size());
+	array<double> var(cnd_col.size());
 
     for (size_t i = 0; i < var.size(); i++)
 	{
-        symbol_table.add_variable(cnd_tars[i], var[i]);
+        symbol_table.add_variable(cnd_col[i], var[i]);
 	}
 
     exprtk::expression<double> expression;
@@ -863,34 +874,79 @@ void gctl::dsv_io::filter(std::string cnd_str, const std::vector<std::string> &c
 	if (!parser.compile(cnd_str, expression))
     {
         throw std::runtime_error("[gctl::dsv_io] Fail to compile the math expression.");
-    }
+    }   
 
-    if (thead == RowHead) // cnd_tars是行头 此时为按列过滤
+    std::vector<std::string> str_line, row_names;
+    std::vector<std::vector<std::string> > str_table; 
+
+    for (size_t i = 1; i <= row_num_; i++)
     {
-        for (size_t i = 1; i <= col_num_; i++)
+        for (size_t j = 0; j < var.size(); j++)
         {
-            for (size_t j = 0; j < var.size(); j++)
+            var[j] = table_[i][idx[j]].value<double>();
+        }
+
+        if (expression.value() > 0.5) // return 1 if matched or 0 if dismatched
+        {
+            if (out_row)
             {
-                var[j] = table_[idx[j]][i].value<double>();
+                str_line.clear();
+
+                for (size_t j = 1; j <= col_num_; j++)
+                {
+                    str_line.push_back(table_[i][j].str_);
+                }
+                
+                str_table.push_back(str_line);
+            }
+            else
+            {
+                str_line.clear();
+                for (size_t j = 0; j < idx.size(); j++)
+                {
+                    str_line.push_back(table_[i][idx[j]].str_);
+                }
+                
+                for (size_t j = 0; j < odx.size(); j++)
+                {
+                    str_line.push_back(table_[i][odx[j]].str_);
+                }
+                
+                str_table.push_back(str_line);
             }
 
-            // return 1 if matched or 0 if dismatched
+            row_names.push_back(table_[i][0].str_);
-            if (expression.value() < 0.5) column_output(i, Disable);
         }
     }
-    else // cnd_tars是列头 此时为按行过滤
-    {
-        for (size_t i = 1; i <= row_num_; i++)
-        {
-            for (size_t j = 0; j < var.size(); j++)
-            {
-                var[j] = table_[i][idx[j]].value<double>();
-            }
 
-            // return 1 if matched or 0 if dismatched
+    out_table.init_table(str_table);
-            if (expression.value() < 0.5) row_output(i, Disable);
+
+    std::vector<std::string> io_col;
+    if (out_row)
+    {
+        get_column_names(io_col);
+        out_table.cell(table_[0][0].str_, 0, 0);
+    }
+    else
+    {
+        for (size_t j = 0; j < idx.size(); j++)
+        {
+            io_col.push_back(cnd_col[j]);
+        }
+
+        for (size_t j = 0; j < odx.size(); j++)
+        {
+            io_col.push_back(out_col[j]);
         }
     }
+
+    out_table.set_column_names(io_col);
+    out_table.set_row_names(row_names, {}, table_[0][0].str_);
+
+    destroy_vector(row_names);
+    destroy_vector(io_col);
+    destroy_vector(str_line);
+    destroy_vector(str_table);
     return;
 }
 
@@ -900,8 +956,9 @@ gctl::geodsv_io::geodsv_io(){}
 
 gctl::geodsv_io::~geodsv_io(){}
 
-gctl::geodsv_io::geodsv_io(std::string filename, std::string file_exten, int t)
+gctl::geodsv_io::geodsv_io(std::string filename, std::string file_exten, table_headtype_e t)
 {
+    file_ = "";
     att_sym_ = '#';
 	tag_sym_ = '!';
 	deli_sym_ = ' ';

lib/io/dsv_io.h

@@ -28,23 +28,19 @@
 #ifndef _GCTL_DSV_IO_H
 #define _GCTL_DSV_IO_H
 
+#include "../gctl_config.h"
+#include "../core.h"
+#include "../utility.h"
+#include "../geometry.h"
 #include "regex.h"
 
-#include "../core/matrix.h"
-#include "../poly/point2c.h"
-#include "../poly/point2p.h"
-#include "../poly/point3c.h"
-#include "../poly/point3s.h"
-#include "file_io.h"
-#include "term_io.h"
-
 #ifdef GCTL_EXPRTK
 #include "exprtk.hpp"
 #endif // GCTL_EXPRTK
 
 namespace gctl
 {
-	enum table_cell_type
+	enum cell_type_e
 	{
 		String,
 		Int,
@@ -54,7 +50,7 @@ namespace gctl
 	struct table_cell
 	{
 		std::string str_; // 单元格的内容 统一保存为字符串
-		table_cell_type type_; // 类型字符串
+		cell_type_e type_; // 类型字符串
 		bool out_ok_; // 是否可输出到文件
 
 		table_cell()
@@ -105,25 +101,15 @@ namespace gctl
 	};
 
 	/**
-	 * @brief 表格的头名称类型
+	 * @brief 表格的头信息类型
+	 * 
 	 */
 	enum table_headtype_e
 	{
-		NoHead = 1, // 0001 无表头
+		NoHead, // 没有表头
-		ColHead = 2, // 0010 有列表头
+		BothHead, // 同时有行与列表头 
-		RowHead = 4, // 0100 有行表头
+		ColumnHead, // 只有列表头
-	};
+		RowHead, // 只有行表头
-
-	/**
-	 * @brief 表格的信息类型
-	 */
-	enum table_infotype_e
-	{
-		TagInfo = 1, // 00001
-		AttInfo = 2, // 00010
-		HeadInfo = 4, // 00100
-		RowInfo = 8, // 01000
-		ColInfo = 16, // 10000
 	};
 
 	/**
@@ -140,6 +126,8 @@ namespace gctl
     class dsv_io
 	{
 	protected:
+		// 文件名
+		std::string file_;
         // 头信息行数 表格行数（不包括表头） 表格列数（不包括表头）
 		int head_num_, row_num_, col_num_;
 		// 注释行起始符 标记行起始符 分割符
@@ -169,13 +157,21 @@ namespace gctl
 		 * @param file_exten 文件扩展名
 		 * @param t 表格是否有行和列名称
 		 */
-		dsv_io(std::string filename, std::string file_exten = ".txt", int t = NoHead);
+		dsv_io(std::string filename, std::string file_exten = ".txt", table_headtype_e t = NoHead);
 
 		/**
-         * @brief 清理（还原）表格
+         * @brief 清理字符串向量对象
+         * 
          */
 		void clear();
 
+		/**
+		 * @brief 返回头信息行数
+		 * 
+		 * @return 行数
+		 */
+		int head_number(){return head_num_;}
+
 		/**
 		 * @brief 返回数据行数
 		 * 
@@ -190,103 +186,89 @@ namespace gctl
 		 */
 		int col_number(){return col_num_;}
 
-		/**
-		 * @brief 返回头信息行数
-		 * 
-		 * @return 行数
-		 */
-		int head_number(){return head_num_;}
-
-		/**
-		 * @brief 设置头信息行数
-		 * 
-		 * @param num 行数 
-		 */
-		void head_number(int num){head_num_ = num;}
-
 		/**
 		 * @brief 返回头信息
 		 * 
 		 * @return 头信息
 		 */
-		const std::vector<std::string> &head_records(){return heads_;}
+		const std::vector<std::string>& get_head_records(){return heads_;}
-
-		/**
-		 * @brief 设置头信息
-		 * 
-		 * @param heads 头信息
-		 */
-		void head_records(const std::vector<std::string> &heads){heads_ = heads; head_num_ = heads_.size();}
-
-		/**
-		 * @brief 设置列分隔符
-		 * 
-		 * @param deli_sym 分隔符
-		 */
-		void delimeter(char deli_sym){deli_sym_ = deli_sym;}
-
-		/**
-		 * @brief 获取列分隔符
-		 *
-		 * @return 分隔符
-		 */
-		char delimeter(){return deli_sym_;}
-
-		/**
-		 * @brief 设置注释行符号
-		 * 
-		 * @param att_sym 注释符号
-		 */
-		void annotation_symbol(char att_sym){att_sym_ = att_sym;}
-
-		/**
-		 * @brief 获取注释行符号
-		 *
-		 * @return 注释符号
-		 */
-		char annotation_symbol(){return att_sym_;}
-
-		/**
-		 * @brief 设置注释
-		 * 
-		 * @param att 注释
-		 */
-		void annotations(const std::vector<std::string> &att){annotates_ = att;}
 
 		/**
 		 * @brief 返回注释行
 		 * 
 		 * @return 注释行
 		 */
-		const std::vector<std::string> &annotations(){return annotates_;}
+		const std::vector<std::string>& get_annotoations(){return annotates_;}
-
-		/**
-		 * @brief 设置标记行符号
-		 * 
-		 * @param tag_sym 标记符号
-		 */
-		void tag_symbol(char tag_sym){tag_sym_ = tag_sym;}
-
-		/**
-		 * @brief 获取标记行符号
-		 *
-		 * @return 标记符号
-		 */
-		char tag_symbol(){return tag_sym_;}
-
-		/**
-		 * @brief 设置标记
-		 * 
-		 * @param tags 标记
-		 */
-		void tags(const std::vector<std::string> &tags){tags_ = tags;}
 
 		/**
 		 * @brief 返回标记行
 		 * 
 		 * @return 标记
 		 */
-		const std::vector<std::string> &tags(){return tags_;}
+		const std::vector<std::string>& get_tags(){return tags_;}
+
+		/**
+		 * @brief 获取行名称
+		 * 
+		 * @param names 名称
+		 */
+		void get_row_names(std::vector<std::string> &names);
+
+		/**
+		 * @brief 获取列名称
+		 * 
+		 * @param names 名称
+		 */
+		void get_column_names(std::vector<std::string> &names);
+
+		/**
+		 * @brief 设置列分隔符
+		 * 
+		 * @param deli_sym 分隔符
+		 */
+		void set_delimeter(char deli_sym){deli_sym_ = deli_sym;}
+
+		/**
+		 * @brief 设置头信息行数
+		 * 
+		 * @param num 行数 
+		 */
+		void set_head_number(char num){head_num_ = num;}
+
+		/**
+		 * @brief 设置注释行符号
+		 * 
+		 * @param att_sym 注释符号
+		 */
+		void set_annotation_symbol(char att_sym){att_sym_ = att_sym;}
+
+		/**
+		 * @brief 设置标记行符号
+		 * 
+		 * @param tag_sym 标记符号
+		 */
+		void set_tag_symbol(char tag_sym){tag_sym_ = tag_sym;}
+
+		/**
+		 * @brief 设置头信息
+		 * 
+		 * @param heads 头信息
+		 */
+		void set_head_records(const std::vector<std::string> &heads){heads_ = heads; head_num_ = heads_.size();}
+
+		/**
+		 * @brief 设置注释
+		 * 
+		 * @param att 注释
+		 */
+		void set_annotoations(const std::vector<std::string> &att){annotates_ = att;}
+
+		/**
+		 * @brief 设置标记
+		 * 
+		 * @param tags 标记
+		 */
+		void set_tags(const std::vector<std::string> &tags){tags_ = tags;}
 
 		/**
 		 * @brief 设置行名称
@@ -295,14 +277,7 @@ namespace gctl
 		 * @param idx 索引数组（若为空则依次设置行名称）
 		 * @param corner_name 表格左上角位置名称（默认为RowNames）
 		 */
-		void row_names(const std::vector<std::string> &names, const std::vector<int> &idx = {}, std::string corner_name = "RowNames");
+		void set_row_names(const std::vector<std::string> &names, const std::vector<int> &idx = {}, std::string corner_name = "RowNames");
-
-		/**
-		 * @brief 获取行名称
-		 * 
-		 * @return 行名称
-		 */
-		std::vector<std::string> row_names();
 
 		/**
 		 * @brief 设置列名称
@@ -310,14 +285,7 @@ namespace gctl
 		 * @param names 名称数组
 		 * @param idx 索引数组（若为空则依次设置行名称）
 		 */
-		void column_names(const std::vector<std::string> &names, const std::vector<int> &idx = {});
+		void set_column_names(const std::vector<std::string> &names, const std::vector<int> &idx = {});
-
-		/**
-		 * @brief 获取列名称
-		 * 
-		 * @return 列名称
-		 */
-		std::vector<std::string> column_names();
 
 		/**
 		 * @brief 设置行类型
@@ -325,7 +293,7 @@ namespace gctl
 		 * @param t 类型名称 String|Float|Int
 		 * @param idx 行索引
 		 */
-		void row_type(table_cell_type t, int idx);
+		void set_row_type(cell_type_e t, int idx);
 
 		/**
 		 * @brief 设置行类型
@@ -333,23 +301,7 @@ namespace gctl
 		 * @param t 类型名称 String|Float|Int
 		 * @param name 行名称
 		 */
-		void row_type(table_cell_type t, std::string name);
+		void set_row_type(cell_type_e t, std::string name);
-
-		/**
-		 * @brief 获取行类型
-		 *
-		 * @param idx 行索引
-		 * @return 行类型
-		 */
-		table_cell_type row_type(int idx);
-
-		/**
-		 * @brief 获取行类型
-		 *
-		 * @param name 行名称
-		 * @return 行类型
-		 */
-		table_cell_type row_type(std::string name);
 
 		/**
 		 * @brief 设置列类型
@@ -357,7 +309,7 @@ namespace gctl
 		 * @param t 类型名称 String|Float|Int
 		 * @param idx 列索引
 		 */
-		void column_type(table_cell_type t, int idx);
+		void set_column_type(cell_type_e t, int idx);
 
 		/**
 		 * @brief 设置列类型
@@ -365,23 +317,7 @@ namespace gctl
 		 * @param t 类型名称 String|Float|Int
 		 * @param name 列名称
 		 */
-		void column_type(table_cell_type t, std::string name);
+		void set_column_type(cell_type_e t, std::string name);
-
-		/**
-		 * @brief 获取列类型
-		 *
-		 * @param idx 列索引
-		 * @return 列类型
-		 */
-		table_cell_type column_type(int idx);
-
-		/**
-		 * @brief 获取列类型
-		 *
-		 * @param name 列名称
-		 * @return 列类型
-		 */
-		table_cell_type column_type(std::string name);
 
 		/**
 		 * @brief 读入文本文件
@@ -390,14 +326,14 @@ namespace gctl
 		 * @param file_exten 文件扩展名
 		 * @param t 表格是否有行和列名称
 		 */
-		void load_text(std::string filename, std::string file_exten = ".txt", int t = NoHead);
+		void load_text(std::string filename, std::string file_exten = ".txt", table_headtype_e t = NoHead);
 
 		/**
 		 * @brief 读入CSV文件
 		 * 
 		 * @param filename 文件名
 		 */
-		void load_csv(std::string filename, int t = ColHead);
+		void load_csv(std::string filename, table_headtype_e t = ColumnHead);
 		
 		/**
 		 * @brief 将内容写入文件
@@ -421,44 +357,23 @@ namespace gctl
 		 * @param col 数据列数
 		 */
 		void init_table(int row, int col);
-		
-		/**
-		 * @brief 导出表格（默认不会导出失效的行和列）
-		 * 
-		 * @param ignore_disabled 忽略失效的行和列
-		 * @return 导出的表格
-		 */
-		dsv_io export_table(bool ignore_disabled = true);
 
 		/**
 		 * @brief 返回表格信息
 		 * 
 		 * @param t 显示表格信息的类型
 		 */
-		void info(int t, std::ostream &os = std::cout);
+		void info(table_headtype_e t = ColumnHead);
-
-		/**
-		 * @brief 显示表格（不会显示失效的行和列）
-		 * 
-		 */
-		void display();
 
 		/**
 		 * @brief 返回名称为name和R<id>和C<id>的行或列的索引
 		 * 
 		 * @param name 名称 可以是具体的名称（如有），或者是R<id>和C<id>
-		 * @param iter_row 搜索行名称（默认为搜索列名称）
+		 * @param iter_row 搜索行名称（默认为搜索列名称），如果name参数为R<id>和C<id>则此参数无效
 		 * @return 索引 返回的索引（大于等于1 小于等于行数或列数）失败则返回-1
 		 */
 		int name_index(std::string name, bool iter_row = false);
 
-		/**
-		 * @brief 设置表格输出（作用于整个表格）
-		 * 
-		 * @param s 设置输出类型
-		 */
-		void table_output(switch_type_e s);
-
 		/**
 		 * @brief 设置列输出。你仍然可以使用这些数据，它们只是不会被输出
 		 * 
@@ -532,27 +447,32 @@ namespace gctl
 		int add_row(std::string name, std::string id_name);
 
 		/**
-		 * @brief 按行或列过滤数据（不符合条件的行与列的输出属性将设置为假）
+		 * @brief 按行过滤并返回符合条件的列数据
+		 * 
+		 * @note  过滤后的表格第一列尾用于匹配正则表达式的列，剩余列尾为筛选后符合条件的列数据。
 		 * 
 		 * @param cnd_str 正则表达式
-		 * @param cnd_tar 用于匹配正则表达式的行或列名称
+		 * @param cnd_col 用于匹配正则表达式的列名称
-		 * @param thead   用于匹配正则表达式的行或列类型 RowHead表示名称为行头则按列过滤，ColHead时表示名称为列头则按行过滤
+		 * @param out_col 输出的列索引列表（列表为空时则会输出所有列），正则表达式为真时即筛选这些行与列上对应的数据
+		 * @param out_table 输出的表格
 		 */
-		void filter(std::string cnd_str, std::string cnd_tar, table_headtype_e thead = RowHead);
+		void filt_column(std::string cnd_str, std::string cnd_col, 
+			const std::vector<std::string> &out_col, dsv_io &out_table);
 
 		/**
-		 * @brief table line operate function pointer
+		 * @brief row operate function pointer
 		 * 
 		 */
-		typedef bool (*linebool_func_t)(const std::vector<table_cell> &table_line);
+		typedef bool (*rowbool_func_t)(const std::vector<table_cell> &table_row);
 
 		/**
-		 * @brief  逐行或列过滤数据（不符合条件的行与列的输出属性将设置为假）
+		 * @brief  按行过滤并返回符合条件的列数据
 		 * 
-		 * @param func  处理行或列类容的布尔函数
+		 * @param func 处理行类容的布尔函数
-		 * @param thead RowHead时表示按行过滤，ColHead时表示按列过滤
+		 * @param out_col 输出的列索引列表（列表为空时则会输出所有列），正则表达式为真时即筛选这些行与列上对应的数据
+		 * @param out_table 输出的表格
 		 */
-		void filter(linebool_func_t func, table_headtype_e thead = RowHead);
+		void filt_column(rowbool_func_t func, const std::vector<std::string> &out_col, dsv_io &out_table);
 
 #ifdef GCTL_EXPRTK
 
@@ -569,16 +489,18 @@ namespace gctl
 		void cal_column(std::string expr_str, const std::vector<std::string> &col_list, int p = 6);
 
 		/**
-		 * @brief 按行或列过滤数据（不符合条件的行与列的输出属性将设置为假）
+		 * @brief 按行过滤并返回符合条件的列数据
 		 * 
 		 * @note  只有单元格类型为float和Int类型的列数据才能用于计算。计算由exprtk库完成，支持的表达式见其说明文档。
 		 * 因为没有使用strtk库的相关内容，所以并不支持对字符串与数字类型的混合条件判断。基于字符串的内容提取请使用其他函数。
 		 * 
 		 * @param cnd_str 条件表达式
-		 * @param cnd_tars 用于计算条件表达式的行或列名称数组
+		 * @param cnd_col 用于条件表达式的列索引列表
-		 * @param thead   用于计算条件表达式的行或列类型 RowHead表示名称为行头则按列过滤，ColHead时表示名称为列头则按行过滤
+		 * @param out_col 输出的列索引列表（列表为空时则会输出所有列），即条件判断为真时即筛选这些行与列上对应的数据
+		 * @param out_table 输出的表格
 		 */
-		void filter(std::string cnd_str, const std::vector<std::string> &cnd_tars, table_headtype_e thead = RowHead);
+		void filt_column(std::string cnd_str, const std::vector<std::string> &cnd_col, 
+			const std::vector<std::string> &out_col, dsv_io &out_table);
 
 #endif // GCTL_EXPRTK
 
@@ -883,7 +805,7 @@ namespace gctl
 		 * @param filename 文件名
 		 * @param file_exten 文件扩展名
 		 */
-		geodsv_io(std::string filename, std::string file_exten = ".txt", int t = NoHead);
+		geodsv_io(std::string filename, std::string file_exten = ".txt", table_headtype_e t = NoHead);
 
 		/**
 		 * @brief 填充二维坐标列
@@ -1045,6 +967,6 @@ namespace gctl
 		 */
 		void get_column_point3ds(array<point3ds> &data, std::string rname, std::string pname, std::string tname);
 	};
-};
+}
 
 #endif //_GCTL_DSV_IO_H

lib/io/gmsh_io.h

@@ -29,21 +29,9 @@
 #define _GCTL_GMSH_IO_H
 
 // library's head files
-#include "../core/matrix.h"
+#include "../core.h"
-#include "../poly/node.h"
+#include "../geometry.h"
-#include "../poly/sphere.h"
+#include "../utility.h"
-#include "../poly/triangle2d.h"
-#include "../poly/triangle2d2o.h"
-#include "../poly/rectangle2d.h"
-#include "../poly/edge2d.h"
-#include "../poly/block.h"
-#include "../poly/prism.h"
-#include "../poly/tesseroid.h"
-#include "../poly/tetrahedron.h"
-#include "../poly/tri_cone.h"
-#include "../poly/triangle.h"
-#include "../poly/edge.h"
-#include "file_io.h"
 
 #ifdef GCTL_EIGEN
 #include "Eigen/Dense"
@@ -113,7 +101,6 @@ namespace gctl
 		template <typename E, typename A> void save_mesh(const array<type_block<E>> &element, const array<vertex<point3dc, A>> &node);
 		template <typename E, typename A> void save_mesh(const array<type_tricone<E>> &element, const array<vertex<point3dc, A>> &node);
 		template <typename E, typename A> void save_mesh(const array<type_prism<E>> &element, const array<vertex<point3dc, A>> &node);
-		template <typename E, typename A> void save_mesh(const array<type_tesseroid<E>> &element, const array<vertex<point3dc, A>> &node);
 
 		template <typename E> void save_mesh(const array<type_tetrahedron<E>> &element);
 
@@ -423,14 +410,6 @@ namespace gctl
 		return;
 	}
 
-	template <typename E, typename A>
-	void gmshio::save_mesh(const array<type_tesseroid<E>> &element, const array<vertex<point3dc, A>> &node)
-	{
-		initialized(Output);
-		save2gmsh(gmsh_out, element, node, out_packed);
-		return;
-	}
-
 	template <typename E>
 	void gmshio::save_mesh(const array<type_tetrahedron<E>> &element)
 	{
@@ -3858,6 +3837,5 @@ namespace gctl
 
 #endif // GCTL_EIGEN
 
-};
+}
-
 #endif //_GCTL_GMSH_IO_H

lib/io/mesh_io.h

@@ -0,0 +1,526 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_MESH_IO_H
+#define _GCTL_MESH_IO_H
+
+// library's head files
+#include "../core.h"
+#include "../geometry.h"
+#include "../utility.h"
+
+#include "triangle_io.h"
+#include "tetgen_io.h"
+#include "gmsh_io.h"
+#include "map"
+
+
+namespace gctl
+{
+    /**
+     * @brief 无效的索引缺省值。
+     */
+    #define DEFAULT_INVALID_TAG -9999
+
+    /**
+     * @brief 网格单元体名称枚举类型。
+     */
+    enum element_type_enum
+    {
+        NotSet,
+        _2NodeLine,
+        _3NodeTriangle,
+        _4NodeQuadrangle,
+        _4NodeTetrahedron,
+        _8NodeHexahedron,
+        _6NodePrism,
+        _5NodePyramid,
+        _3NodeSecondOrderLine,
+        _6NdoeSecondOrderLine,
+        _9NodeSecondOrderQuadrangle,
+        _10NodeSecondOrderTetrahedron,
+        _27NodeSecondOrderHexahedron,
+        _18NodeSecondOrderPrism,
+        _14NodeSecondOrderPyramid,
+        _1NodePoint,
+        _8NodeSecondOrderQuadrangle,
+        _20NdoeSecondOrderHexahedron,
+        _15NodeSecondOrderPrism,
+        _13NodeSecondOrderPyramid,
+        _9NodeThirdOrderIncompleteTriangle,
+        _10NdoeThirdOrderTriangle,
+        _12NodeFourthOrderIncompleteTriangle,
+        _15NodeFourthOrderTriangle,
+        _15NodeFifthOrderCompleteTriangle,
+        _21NodeFifthOrderCompleteTriangle,
+        _4NodeThirdOrderEdge,
+        _5NodeFourthOrderEdge,
+        _6NodeFifthOrderEdge,
+        _20NodeThirdOrderTetrahedron,
+        _35NodeFourthOrderTetrahedron,
+        _56NodeFifithOrderTetrahedron,
+        _64NodeThirdOrderHexahedron,
+        _125NodeFourthOrderHexahedron,
+    };
+
+    /**
+     * @brief 网格单元体标签类型枚举类型
+     * 
+     */
+    enum element_tag_enum
+    {
+        PhysicalTag, // 元素的物理分组标签
+        GeometryTag, // 元素的几何分组标签
+        PartitionTag, // 元素的剖分分组标签
+        NodeTag, // 顶点的标签（仅用于输出顶点标签数据）
+    };
+
+    /**
+     * @brief 网格单元体结构体
+     * 
+     */
+    struct meshio_element
+    {
+        bool enabled; // 单元体是否有效
+        int id; // 单元体编号
+        element_type_enum type; // 单元体类型
+        array<vertex3dc*> vert_ptrs; // 顶点指针数组
+        array<meshio_element*> neigh_ptrs; // 相邻单元体指针数组
+
+        meshio_element();
+    };
+
+    /**
+     * @brief 网格数据结构体
+     * 
+     */
+    struct meshio_data
+    {
+        bool enabled; // 数据体是否有效
+        mesh_data_type_e d_type; // 数据类型
+        array<std::string> str_tag; // 字符串类型的标签（默认为一个，即为数据的名称）
+        array<double> real_tag; // 实数类型的标签（默认为一个，等于0.0）
+        array<int> int_tag; // 整数类型的标签（最少三个，最后一个为数据的长度）
+        array<void*> tar_ptrs; // 数据连接的对象指针数组 具体类型为vertex3dc*或meshio_element*
+        array<double> val; // 数据值（目前仅支持标量数据，后续再添加对矢量数据的支持）
+
+        meshio_data();
+
+        /**
+         * @brief 清空数组并重置变量。
+         * 
+         */
+        void clear();
+
+        /**
+         * @brief 检查数据体是否合规
+         * 
+         */
+        bool pass_check();
+    };
+
+    /**
+     * @brief 网格单元体分组结构体。
+     * 
+     */
+    struct meshio_element_group
+    {
+        bool enabled; // 组是否有效
+        element_type_enum type; // 组内单元体类型
+        int phys_group; // 物理分组标签
+        int geom_group; // 几何分组标签
+        int part_group; // 剖分分组标签
+        std::string name; // 组名
+        std::vector<meshio_element*> elem_ptrs; // 组内单元体指针数组
+
+        meshio_element_group();
+
+        /**
+         * @brief 将组内所有单元体设置为有效状态。
+         * 
+         */
+        void enable_elements();
+
+        /**
+         * @brief 将组内所有单元体设置为无效状态。
+         * 
+         */
+        void disable_elements();
+    };
+    
+    /**
+     * @brief 网格读写类，这个类实现了多种数据格式的网格文件的读写操作。并具备简单的单元体操作功能。
+     * 
+     */
+    class mesh_io
+    {
+    public:
+        mesh_io();
+        virtual ~mesh_io();
+
+        /**
+         * @brief 重置（清空）网格数据。
+         * 
+         */
+        void reset();
+
+        /**
+         * @brief 输出网格数据信息至指定流。
+         * 
+         * @param ss 指定流（默认为clog）
+         */
+        void info(std::ostream &ss = std::clog);
+
+        /**
+         * @brief 编辑网格数据状态。
+         * 
+         * @param swt 使能类型（Enable或Disable）
+         * @param dataname 数据名称（缺省值为null，表示对所有数据进行操作）。
+         */
+        void edit_data(switch_type_e swt, std::string dataname = "null");
+
+        /**
+         * @brief 按单元体类型编辑网格单元体组。
+         * 
+         * @param swt 使能类型（Enable或Disable）
+         * @param e_type 单元体类型（缺省值值NotSet，表示对所有单元体组进行操作）。
+         */
+        void edit_group(switch_type_e swt, element_type_enum e_type = NotSet);
+
+        /**
+         * @brief 按单元体组名称编辑网格单元体组。
+         * 
+         * @param swt 使能类型（Enable或Disable）。
+         * @param grp_name 单元体组名称。
+         */
+        void edit_group(switch_type_e swt, std::string grp_name);
+
+        /**
+         * @brief 按单元体组标签编辑网格单元体组。
+         * 
+         * @param swt 使能类型（Enable或Disable）。
+         * @param tag_type 标签类型（PhysicalTag，GeometryTag或者PartitionTag）。
+         * @param tag 标签值。
+         */
+        void edit_group(switch_type_e swt, element_tag_enum tag_type, int tag);
+
+        /**
+         * @brief 按单元体组标签编辑网格单元体组的名称。
+         * 
+         * @param anchor_type 搜索的标签类型（PhysicalTag，GeometryTag或者PartitionTag）。
+         * @param anchor_group 搜索的标签值。
+         * @param new_name 单元体组的新名称。
+         */
+        void edit_group(element_tag_enum anchor_type, int anchor_group, std::string new_name);
+
+        /**
+         * @brief 按单元体组标签搜索并编辑网格单元体组的标签。
+         * 
+         * @param anchor_type 搜索的标签类型（PhysicalTag，GeometryTag或者PartitionTag）。
+         * @param anchor_group 搜索的标签值
+         * @param tar_type 更改的标签类型（PhysicalTag，GeometryTag或者PartitionTag）。
+         * @param tar_group 更改的标签值
+         */
+        void edit_group(element_tag_enum anchor_type, int anchor_group, element_tag_enum tar_type, int tar_group);
+
+        /**
+         * @brief 返回顶点标签数组的引用。
+         * 
+         * @return 整型数组的引用。
+         */
+        const array<int> &get_node_tag();
+
+        /**
+         * @brief 返回指定类型与名称的标签值
+         * 
+         * @param anchor_type 查找的标签类型（PhysicalTag，GeometryTag或者PartitionTag）。
+         * @param anchor_name 查找的元素组名称。
+         * @return 标签值
+         */
+        int get_tag(element_tag_enum anchor_type, std::string anchor_name);
+
+        /**
+         * @brief 返回所有顶点数组的引用。
+         * 
+         * @return 顶点数组的引用。
+         */
+        const array<vertex3dc> &get_nodes();
+
+        /**
+         * @brief 返回所有单元体数组的引用。
+         * 
+         * @return 单元体数组的引用。
+         */
+        const array<meshio_element> &get_elems();
+
+        /**
+         * @brief 返回对应类型单元体的数量（注意只会统计有效的单元体组）。
+         * 
+         * @param e_type 单元体类型（缺省为NotSet，返回所有单元体类型的总和）。
+         * @return 整型大小。
+         */
+        size_t element_size(element_type_enum e_type = NotSet);
+
+        /**
+         * @brief 返回对应标签类型与标签值的单元体数量（注意只会统计有效的单元体组）。
+         * 
+         * @param tag_type 标签类型。
+         * @param tag 标签值。
+         * @return 整型大小。
+         */
+        size_t element_size(element_tag_enum tag_type, int tag);
+
+        /**
+         * @brief 返回对应名称的单元体数量（注意只会统计有效的单元体组）。
+         * 
+         * @param phys_name 单元体组名称
+         * @return 整型大小。
+         */
+        size_t element_size(std::string phys_name);
+
+        /**
+         * @brief 将对应标签类型转换为网格数据类型
+         * 
+         * @param tag_type 标签类型。
+         */
+        void convert_tags_to_data(element_tag_enum tag_type);
+
+        /**
+         * @brief 输出对应名称的单元体组到三角形数组。
+         * 
+         * @param tris 输出的三角形数组的引用。
+         * @param phys_name 单元体组的名称。
+         */
+        void export_elements_to(array<triangle> &tris, std::string phys_name = "All");
+
+        /**
+         * @brief 输出对应标签类型与标签值的单元体组到三角形数组。
+         * 
+         * @param tris 输出的三角形数组的引用。
+         * @param tag_type 标签类型。
+         * @param tag 标签值。
+         */
+        void export_elements_to(array<triangle> &tris, element_tag_enum tag_type, int tag);
+
+        /**
+         * @brief 输出对应名称的单元体组到四面体数组。
+         * 
+         * @param tets 输出的四面体数组的引用。
+         * @param phys_name 单元体组的名称。
+         */
+        void export_elements_to(array<tetrahedron> &tets, std::string phys_name = "All");
+
+        /**
+         * @brief 输出对应标签类型与标签值的单元体组到四面体数组。
+         * 
+         * @param tris 输出的四面体数组的引用。
+         * @param tag_type 标签类型。
+         * @param tag 标签值。
+         */
+        void export_elements_to(array<tetrahedron> &tets, element_tag_enum tag_type, int tag);
+
+        /**
+         * @brief 获取gmsh格式分组表
+         * 
+         * @param g_groups gmsh格式表
+         */
+        void get_gmsh_physical_groups(std::vector<gmsh_physical_group> &g_groups);
+
+        /**
+         * @brief 检查是否存在名为name的数据
+         * 
+         * @param name 数据名称
+         * @param type 数据类型
+         * 
+         * @return 存在则返回数据索引，不存在则返回-1。
+         */
+        int if_saved_data(std::string name, mesh_data_type_e type);
+
+        /**
+         * @brief 获取数据对象的引用
+         * 
+         * @param name 数据名称
+         * @param type 数据类型
+         * @return 数据引用 
+         */
+        meshio_data &get_data(std::string name, mesh_data_type_e type);
+
+        /**
+         * @brief 获取数据对象的指针
+         * 
+         * @param name 数据名称
+         * @param type 数据类型
+         * @return 数据指针
+         */
+        meshio_data *get_data_ptr(std::string name, mesh_data_type_e type);
+
+        /**
+         * @brief 添加一个顶点数据对象。数据将依次添加到所有顶点位置。
+         * 
+         * @note 若对应名称的数据已经存在则会覆盖
+         * 
+         * @param data 输入的数据数组，长度与网格所有顶点数据相同。
+         * @param name 新建的数据名称。
+         */
+        void add_node_data(std::string name, const array<double> &data);
+
+        /**
+         * @brief 添加一个顶点数据对象。数据将依次添加到布尔为真的顶点位置。
+         * 
+         * @note 若对应名称的数据已经存在则会覆盖
+         * 
+         * @param data 输入的数据数组，长度与网格所有顶点数据相同。
+         * @param boolen 输入的布尔，只有为真元素位置的顶点数据将被保存。
+         * @param name 新建的数据名称。
+         */
+        void add_node_data(std::string name, const array<double> &data, const array<bool> &boolen);
+
+        /**
+         * @brief 按单元体类型筛选创建一个单元体数据对象。数据将依次添加到所选元素位置。
+         * 
+         * @note 若对应名称的数据已经存在则会覆盖
+         * 
+         * @param data 输入的数据数组。
+         * @param name 新建数据名称。
+         * @param e_type 新建数据的单元体类型（缺省值为NotSet，表示选择所有有效的单元体）。
+         */
+        void add_element_data(std::string name, const array<double> &data, element_type_enum e_type = NotSet);
+
+        /**
+         * @brief 按单元体标签值筛选创建一个单元体数据对象。
+         * 
+         * @note 若对应名称的数据已经存在则会追加
+         * 
+         * @param data 输入的数据数组。
+         * @param name 新建数据名称。
+         * @param tag_type 标签类型。
+         * @param tag 标签值。
+         */
+        void add_element_data(std::string name, const array<double> &data, element_tag_enum tag_type, int tag);
+
+        /**
+         * @brief 按单元体组的名称筛选创建一个单元体数据对象。
+         * 
+         * @note 若对应名称的数据已经存在则会追加
+         * 
+         * @param data 输入的数据数组。
+         * @param name 新建数据名称。
+         * @param phys_name 单元体组的名称。
+         */
+        void add_element_data(std::string name, std::string phys_name, const array<double> &data);
+
+        /**
+         * @brief 按单元体组的名称筛选创建一个单元体数据对象。
+         * 
+         * @note 若对应名称的数据已经存在则会追加
+         * 
+         * @param phys_val 数据初始值
+         * @param name 新建数据名称。
+         * @param phys_name 单元体组的名称。
+         */
+        void add_element_data(std::string name, std::string phys_name, double phys_val);
+
+        /**
+         * @brief 读入triangle软件输出的网格剖分文件。
+         * 
+         * @param filename 文件名（.node和.ele文件必须在同一路径下，.neigh文件不是必须的，文件名不包含后缀名）。
+         * @param is_packed 输入文件的索引值是否从0开始。
+         */
+        void read_triangle_ascii(std::string filename, index_packed_e is_packed = Packed);
+
+        /**
+         * @brief 读入tetgen软件输出的网格剖分文件。
+         * 
+         * @param filename 文件名（.node和.ele文件必须在同一路径下，.neigh文件不是必须的，文件名不包含后缀名）。
+         * @param is_packed 输入文件的索引值是否从0开始。
+         */
+        void read_tetgen_ascii(std::string filename, index_packed_e is_packed = Packed);
+
+        /**
+         * @brief 读入Gmsh软件输出的网格剖分文件（只支持v2.2的ASCII文件）。
+         * 
+         * @param filename 文件名
+         * @param is_packed 输入文件的索引值是否从0开始。
+         */
+        void read_gmsh_v2_ascii(std::string filename, index_packed_e is_packed = NotPacked);
+
+        /**
+         * @brief 保存Gmsh软件格式的网格剖分文件（只支持v2.2的ASCII文件）。
+         * 
+         * @param filename 文件名
+         * @param is_packed 输入文件的索引值是否从0开始。
+         */
+        void save_gmsh_v2_ascii(std::string filename, index_packed_e is_packed = NotPacked);
+
+        /**
+         * @brief 保存Paraview软件格式的网格剖分文件
+         * 
+         * @param filename 文件名
+         */
+        void save_vtk_legacy_ascii(std::string filename);
+
+        /**
+         * @brief 导出数据到点云文件
+         * 
+         * @param filename 输出文件名
+         * @param dataname 数据名称
+         * @param out_coor 数据的坐标系（Cartesian或Spherical）
+         * @param refr 参考椭球的短半径（out_coor为Cartesian时无效）
+         * @param refR 参考椭球的长半径（out_coor为Cartesian时无效）
+         */
+        void save_data_to_xyz(std::string filename, std::string dataname = "null", coordinate_system_e out_coor = Cartesian, double refr = GCTL_Earth_Radius, double refR = GCTL_Earth_Radius);
+    
+    protected:
+        bool initialized_; // 类型是否已经初始化完成
+        
+        // 有效的顶点、单元体和单元体组的数量
+        size_t valid_node_size_, valid_elem_size_, valid_group_size_;
+        array<vertex3dc> nodes_; // 网格顶点 当顶点索引为无效值时将不会被输出
+        array<meshio_element> elems_; // 网格元素
+        std::vector<meshio_data> datas_; // 网格数据
+        std::vector<meshio_element_group> groups_; // 网格单元体组
+        array<int> nodes_tag_; // 顶点标签
+        
+        element_type_enum elem_gmshcode2type_[94]; // gmsh的单元体类型数组 数组索引为gmsh的单元体类型码值
+        element_type_enum elem_vtkcode2type_[14]; // vtk的单元体类型数组 数组索引为vtk的单元体类型码值
+        std::map<element_type_enum, int> elem_type2gmshcode_; // 单元体类型到gmsh类型码值的映射
+        std::map<element_type_enum, int> elem_type2vtkcode_; // 单元体类型到vtk类型码值的映射
+        std::map<element_type_enum, int> elem_type2size_; // 单元体类型到单元体顶点数量的映射
+        std::map<element_type_enum, std::string> elem_type2name_; // 单元体类型到单元体名称的映射
+
+        std::string elem_name(element_type_enum e_type); // 获取单元体名称字符串
+        int elem_gmsh_code(element_type_enum e_type); // 获取单元体gmsh类型码值
+        int elem_vtk_code(element_type_enum e_type); // 获取单元体vtk类型码值
+        int elem_size(element_type_enum e_type); // 获取单元体顶点数量
+        element_type_enum elem_gmsh_type(int code); // 获取对应gmsh类型码的单元体类型
+        element_type_enum elem_vtk_type(int code); // 获取对应vtk类型码的单元体类型
+        void update_indexing(); // 更新索引（对网格的元素进行操作后需调用）
+        void sort_groups(); // 对单元体组进行梳理（对网格的元素进行操作后需调用）
+    };
+}
+
+#endif // _GCTL_MESH_IO_H

lib/io/native_io.h

@@ -28,8 +28,8 @@
 #ifndef _GCTL_NATIVE_IO_H
 #define _GCTL_NATIVE_IO_H
 
-#include "../core/spmat.h"
+#include "../core.h"
-#include "file_io.h"
+#include "../utility.h"
 
 namespace gctl
 {
@@ -428,6 +428,6 @@ namespace gctl
         outfile.close();
         return;
     }
-};
+}
 
 #endif // _GCTL_NATIVE_IO_H

lib/io/netcdf_io.h

@@ -28,7 +28,10 @@
 #ifndef _GCTL_NETCDF_IO_H
 #define _GCTL_NETCDF_IO_H
 
-#include "../core/matrix.h"
+// library's head file
+#include "../gctl_config.h"
+#include "../core.h"
+#include "../utility.h"
 
 #ifdef GCTL_NETCDF
 
@@ -771,7 +774,7 @@ namespace gctl
 
 		return;
 	}
-};
+}
 
 #endif // GCTL_NETCDF
 

lib/io/off_io.h

@@ -28,13 +28,15 @@
 #ifndef _GCTL_OFF_IO_H
 #define _GCTL_OFF_IO_H
 
-#include "../poly/triangle.h"
+// library's head file
-#include "file_io.h"
+#include "../core.h"
+#include "../utility.h"
+#include "../geometry.h"
 
 namespace gctl
 {
     template <typename E, typename A>
-    void read_off_acsii(std::string filename, array<vertex<point3dc, E>> &nodes, 
+    void read_Geomview_off(std::string filename, array<vertex<point3dc, E>> &nodes, 
         array<type_triangle<A>> &facets)
     {
         std::ifstream infile;
@@ -118,7 +120,7 @@ namespace gctl
     }
 
     template <typename E, typename A>
-    void save_off_ascii(std::string filename, const array<vertex<point3dc, E>> &nodes, 
+    void save_Geomview_off(std::string filename, const array<vertex<point3dc, E>> &nodes, 
         const array<type_triangle<A>> &facets)
     {
         std::ofstream outfile;
@@ -139,6 +141,6 @@ namespace gctl
         outfile.close();
         return;
     }
-};
+}
 
 #endif // _GCTL_OFF_IO_H

lib/io/ply_io.h

@@ -1,214 +0,0 @@
-/********************************************************
- *  ██████╗  ██████╗████████╗██╗
- * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗██║        ██║   ██║
- * ██║   ██║██║        ██║   ██║
- * ╚██████╔╝╚██████╗   ██║   ███████╗
- *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
- * Geophysical Computational Tools & Library (GCTL)
- *
- * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
- *
- * GCTL is distributed under a dual licensing scheme. You can redistribute 
- * it and/or modify it under the terms of the GNU Lesser General Public 
- * License as published by the Free Software Foundation, either version 2 
- * of the License, or (at your option) any later version. You should have 
- * received a copy of the GNU Lesser General Public License along with this 
- * program. If not, see <http://www.gnu.org/licenses/>.
- * 
- * If the terms and conditions of the LGPL v.2. would prevent you from using 
- * the GCTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GCTL's original author. As a rule, 
- * licenses are provided "as-is", unlimited in time for a one time fee. Please 
- * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
- * to include some description of your company and the realm of its activities. 
- * Also add information on how to contact you by electronic and paper mail.
- ******************************************************/
-
-#ifndef _GCTL_PLY_IO_H
-#define _GCTL_PLY_IO_H
-
-#include "../poly/triangle.h"
-#include "file_io.h"
-
-namespace gctl
-{
-    template <typename E, typename A>
-    void read_ply_ascii(std::string filename, array<vertex<point3dc, E>> &nodes, 
-        array<type_triangle<A>> &facets)
-    {
-        std::ifstream infile;
-        open_infile(infile, filename, ".ply");
-
-        int n, f;
-        std::string line;
-        std::stringstream ss;
-        while (std::getline(infile, line))
-        {
-            if (sscanf(line.c_str(), "element vertex %d", &n)==1) nodes.resize(n);
-            if (sscanf(line.c_str(), "element face %d", &f)==1) facets.resize(f);
-            if (line == "end_header")
-            {
-                if (nodes.empty() || facets.empty())
-                {
-                    infile.close();
-                    throw std::runtime_error("[gctl::read_ply_ascii] Invalid node or facet sizes.");
-                }
-
-                for (size_t i = 0; i < nodes.size(); ++i)
-                {
-                    std::getline(infile, line);
-                    sscanf(line.c_str(), "%lf %lf %lf", &nodes[i].x, &nodes[i].y, &nodes[i].z);
-                    nodes[i].id = i;
-                }
-
-                for (size_t i = 0; i < facets.size(); ++i)
-                {
-                    std::getline(infile, line);
-                    int d, v1, v2, v3;
-                    sscanf(line.c_str(), "%d %d %d %d", &d, &v1, &v2, &v3);
-                    facets[i].vert[0] = &nodes[v1];
-                    facets[i].vert[1] = &nodes[v2];
-                    facets[i].vert[2] = &nodes[v3];
-                    facets[i].id = i;
-                }
-                break;
-            }
-        }
-
-        infile.close();
-        return;
-    }
-
-    template <typename E, typename A>
-    void save_ply_ascii(std::string filename, const array<vertex<point3dc, E>> &nodes, 
-        const array<type_triangle<A>> &facets)
-    {
-        time_t now = time(0);
-        char* dt = ctime(&now);
-
-        std::ofstream outfile;
-        open_outfile(outfile, filename, ".ply");
-
-        outfile << "ply" << std::endl;
-        outfile << "format ascii 1.0" << std::endl;
-        outfile << "comment Created by GCTL at " << dt;
-        outfile << "element vertex " << nodes.size() << std::endl;
-        outfile << "property float x" << std::endl;
-        outfile << "property float y" << std::endl;
-        outfile << "property float z" << std::endl;
-        outfile << "element face " << facets.size() << std::endl;
-        outfile << "property list uchar int vertex_indices" << std::endl;
-        outfile << "end_header" << std::endl;
-
-        for (size_t i = 0; i < nodes.size(); ++i)
-        {
-            outfile << nodes[i].x << " " << nodes[i].y << " " << nodes[i].z << std::endl;
-        }
-
-        for (size_t i = 0; i < facets.size(); ++i)
-        {
-            outfile << "3 " << facets[i].vert[0]->id << " " << facets[i].vert[1]->id << " " << facets[i].vert[2]->id << std::endl;
-        }
-        outfile.close();
-        return;
-    }
-
-    template <typename E, typename A>
-    void read_ply_binary(std::string filename, array<vertex<point3dc, E>> &nodes, 
-        array<type_triangle<A>> &facets)
-    {
-        std::ifstream infile;
-        open_infile(infile, filename, ".ply", std::ios::binary);
-
-        int n, f;
-        std::string line;
-        std::stringstream ss;
-        while (std::getline(infile, line))
-        {
-            if (sscanf(line.c_str(), "element vertex %d", &n)==1) nodes.resize(n);
-            if (sscanf(line.c_str(), "element face %d", &f)==1) facets.resize(f);
-            if (line == "end_header")
-            {
-                if (nodes.empty() || facets.empty())
-                {
-                    infile.close();
-                    throw std::runtime_error("[gctl::read_ply_ascii] Invalid node or facet sizes.");
-                }
-
-                float x, y, z;
-                for (size_t i = 0; i < nodes.size(); ++i)
-                {
-                    infile.read((char*)&x, sizeof(float));
-                    infile.read((char*)&y, sizeof(float));
-                    infile.read((char*)&z, sizeof(float));
-                    nodes[i].x = x;
-                    nodes[i].y = y;
-                    nodes[i].z = z;
-                    nodes[i].id = i;
-                }
-
-                unsigned char num_vertices;
-                int v1, v2, v3;
-                for (size_t i = 0; i < facets.size(); ++i)
-                {
-                    infile.read((char*)&num_vertices, sizeof(unsigned char));
-                    infile.read((char*)&v1, sizeof(int));
-                    infile.read((char*)&v2, sizeof(int));
-                    infile.read((char*)&v3, sizeof(int));
-                    facets[i].vert[0] = &nodes[v1];
-                    facets[i].vert[1] = &nodes[v2];
-                    facets[i].vert[2] = &nodes[v3];
-                    facets[i].id = i;
-                }
-                break;
-            }
-        }
-
-        infile.close();
-        return;
-    }
-
-    template <typename E, typename A>
-    void save_ply_binary(std::string filename, const array<vertex<point3dc, E>> &nodes, 
-        const array<type_triangle<A>> &facets)
-    {
-        std::ofstream outfile;
-        open_outfile(outfile, filename, ".ply", std::ios::binary);
-
-        outfile << "ply\n"
-                << "format binary_little_endian 1.0\n"
-                << "comment Created by GCTL\n"
-                << "element vertex " << nodes.size() << "\n"
-                << "property float x\n"
-                << "property float y\n"
-                << "property float z\n"
-                << "element face " << facets.size() << "\n"
-                << "property list uchar int vertex_indices\n"
-                << "end_header\n";
-
-        float x, y, z;
-        for (size_t i = 0; i < nodes.size(); ++i)
-        {
-            x = nodes[i].x;
-            y = nodes[i].y;
-            z = nodes[i].z;
-            outfile.write((char*)&x, sizeof(float));
-            outfile.write((char*)&y, sizeof(float));
-            outfile.write((char*)&z, sizeof(float));
-        }
-
-        unsigned char num_vertices = 3;
-        for (size_t i = 0; i < facets.size(); ++i)
-        {
-            outfile.write((char*)&num_vertices, sizeof(unsigned char));
-            outfile.write((char*)&facets[i].vert[0]->id, sizeof(int));
-            outfile.write((char*)&facets[i].vert[1]->id, sizeof(int));
-            outfile.write((char*)&facets[i].vert[2]->id, sizeof(int));
-        }
-        outfile.close();
-        return;
-    }
-};
-
-#endif // _GCTL_PLY_IO_H