Compare commits

...

55 Commits
main ... dev_yi

Author SHA1 Message Date
8ae228a0d3 tmp 2025-04-24 12:57:06 +08:00
d935a906e0 tmp 2025-04-23 20:59:46 +08:00
c647ae37e9 move mesh_io to mesh 2025-04-23 17:05:28 +08:00
f18ad0c6cf update cmakelists 2025-04-23 13:53:46 +08:00
4030726e33 tmp 2025-04-23 13:43:40 +08:00
11d29c3dd4 tmp 2025-04-23 12:39:44 +08:00
0e1f69c3d3 bug fixed in sparray 2025-04-08 16:23:40 +08:00
6e944f182a tmp 2025-04-04 19:53:32 +08:00
1447066ed4 tmp 2025-03-06 14:21:04 +08:00
45d38f5f50 tmp 2025-03-06 11:00:59 +08:00
17f1b68087 tmp 2025-03-05 22:37:40 +08:00
4e80b85474 tmp 2025-03-05 15:42:47 +08:00
5528c948e1 tmp 2025-03-05 09:19:16 +08:00
3ca661f9e6 tmp 2025-03-04 22:24:20 +08:00
02deea8162 tmp 2025-03-04 12:14:10 +08:00
71068038e8 update dsviewer 2025-03-03 09:15:25 +08:00
8e46c178be tmp 2025-02-27 09:53:00 +08:00
cf1b1ab4e4 update dsviewer 2025-02-27 09:13:23 +08:00
dfce004048 tmp 2025-02-26 21:37:39 +08:00
a43fd84986 tmp 2025-02-24 21:45:40 +08:00
a1bf953d4b tmp 2025-02-24 12:09:53 +08:00
1ddfdf227f tmp 2025-02-23 10:35:11 +08:00
ef6ab8caa7 mesh_io bug fixed 2025-02-22 19:54:06 +08:00
3302548788 tmp 2025-02-22 11:40:45 +08:00
99a2de65d6 update matrix 2025-02-21 09:51:44 +08:00
db71a9a47f tmp 2025-02-20 12:50:16 +08:00
0f4a1090d9 tmp 2025-02-19 10:22:01 +08:00
5a00846fa8 tmp 2025-02-18 22:02:35 +08:00
6c7ec61fca update mesh_io 2025-02-17 11:50:14 +08:00
91830b3fdd tmp 2025-02-16 20:02:31 +08:00
fb4ffd2acd tmp 2025-02-15 21:50:19 +08:00
10c884b155 tmp 2025-02-14 09:45:57 +08:00
908fa45d1b update dsv_io 2025-02-13 21:53:33 +08:00
140848e882 tmp 2025-02-13 16:41:14 +08:00
014d12cec4 update dsv_io 2025-02-13 10:47:37 +08:00
75ad8a3d27 update dsv_io 2025-02-13 10:44:07 +08:00
12e15b869b add type to dsv_io 2025-02-13 10:31:15 +08:00
fe46cc7aa0 update dsv_io 2025-02-13 10:10:19 +08:00
2802086a76 update dsv_io 2025-02-12 13:19:22 +08:00
fc34ed3e2c tmp 2025-02-11 22:20:51 +08:00
b184ac95fb ymp 2025-02-10 21:25:13 +08:00
6569b9a1b9 tmp 2025-02-09 19:08:36 +08:00
86e0dcd50d tmp 2025-02-09 13:12:55 +08:00
428442629d tmp 2025-02-08 21:04:59 +08:00
0212883ec7 tmp 2025-02-07 11:56:53 +08:00
c413814a8b tmp 2025-02-06 21:17:24 +08:00
51bca90178 add functions 2025-02-05 14:25:59 +08:00
2215a17d5f tmp 2025-02-05 14:11:40 +08:00
b011a82f64 tmp 2025-02-04 11:58:22 +08:00
00270a7b0b tmp 2025-02-04 11:56:40 +08:00
e44f34958e tmp 2025-02-04 11:54:29 +08:00
4e7a7f3791 tmp 2025-02-03 22:30:56 +08:00
2ce5da95b0 tmp 2025-02-03 15:41:30 +08:00
8b294f61d6 update 2025-02-03 12:28:45 +08:00
ebab9756c7 tmp 2025-02-03 12:27:35 +08:00
142 changed files with 5143 additions and 5462 deletions

View File

@ -1,9 +1,15 @@
cmake_minimum_required(VERSION 3.15.2)
#
project(GCTL VERSION 1.0)
project(GCTL VERSION 2.0)
#
include(CMakePackageConfigHelpers)
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
# ExprTKmacOS 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)
@ -11,6 +17,7 @@ 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)
#
@ -26,17 +33,26 @@ 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_INC_DIR})
include_directories(${FFTW3_INCLUDE_DIRS})
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)
@ -54,6 +70,8 @@ if(GCTL_OPENBLAS)
endif()
if(GCTL_NETCDF)
find_package(HDF5)
if(NOT netCDF_FOUND)
find_package(netCDF REQUIRED)
include_directories(${netCDF_INCLUDE_DIR})

View File

@ -16,6 +16,7 @@ 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@)
@ -25,12 +26,15 @@ 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})
@ -45,10 +49,17 @@ endif()
if(@PROJECT_NAME@_FFTW3)
if(NOT FFTW3_FOUND)
find_package(FFTW3 REQUIRED)
include_directories(${FFTW3_INC_DIR})
include_directories(${FFTW3_INCLUDE_DIRS})
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)

View File

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

View File

@ -23,10 +23,17 @@ add_example(windowfunc_ex OFF)
add_example(legendre_ex OFF)
add_example(refellipsoid_ex OFF)
add_example(kde_ex OFF)
add_example(meshio_ex ON)
add_example(meshio_ex OFF)
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(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)

47
example/cliplot_ex.cpp Normal file
View File

@ -0,0 +1,47 @@
/********************************************************
*
*
*
*
*
*
* 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;
}

View File

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

65
example/gmt_ex.cpp Normal file
View File

@ -0,0 +1,65 @@
/********************************************************
*
*
*
*
*
*
* 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;
}

97
example/gnuplot_ex.cpp Normal file
View File

@ -0,0 +1,97 @@
/********************************************************
*
*
*
*
*
*
* 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;
}

View File

@ -1,87 +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
{
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);
}

44
example/ply_io_ex.cpp Normal file
View File

@ -0,0 +1,44 @@
/********************************************************
*
*
*
*
*
*
* 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);
}

View File

@ -25,20 +25,48 @@
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#include "../lib/geodesy.h"
#include "../lib/geometry.h"
using namespace gctl;
int main(int argc, char const *argv[]) try
{
refellipsoid relli(WGS84);
refellipsoid relli2(1000, 900, false);
refellipsoid relli3(1000, 10, true);
refellipsoid ellip;
ellip.set(WGS84);
/*
point3ds ps;
ps.set_io_precision(17);
ps.lon = -110;
ellip.geodetic2spherical(-40, 0, ps.lat, ps.rad);
std::cout << "Geocentric: " << ps << std::endl;
for (size_t i = 0; i <= 90; i++)
{
std::cout << i << " " << relli.radis_at(1.0*i) << " " << relli2.radis_at(1.0*i) << " " << relli3.radis_at(1.0*i) << "\n";
}
point3dc pc = ps.s2c();
pc.set_io_precision(17);
std::cout << "XYZ: " << pc << std::endl;
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)
{

View File

@ -7,7 +7,7 @@
*
* Geophysical Computational Tools & Library (GCTL)
*
* Copyright (c) 2023 Yi Zhang (yizhang-geo@zju.edu.cn)
* 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
@ -25,16 +25,19 @@
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#ifndef _GCTL_MATHS_H
#define _GCTL_MATHS_H
#include "gctl/core.h"
#include "maths/mathfunc_t.h"
#include "maths/mathfunc.h"
#include "maths/mathfunc_ext.h"
#include "maths/shapefunc.h"
#include "maths/legendre.h"
#include "maths/linear_algebra.h"
#include "maths/gaussfunc.h"
#include "maths/fft.h"
int main(int argc, char const *argv[])
{
gctl::sparray2d<double> M(10, 10, 0.0);
#endif // _GCTL_MATHS_H
M.at(4)->set(4, 5);
M.at(2)->set(3, 8);
M.at(6)->set(8, 2);
for (int i = 0; i < M.size(); i++)
{
M.at(i)->show_full();
}
return 0;
}

114
example/sparray_ex.cpp Normal file
View File

@ -0,0 +1,114 @@
/********************************************************
*
*
*
*
*
*
* 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;
}

View File

@ -1,4 +1,4 @@
/********************************************************
/********************************************************
*
*
*
@ -25,24 +25,21 @@
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#ifndef _GCTL_ALGORITHM_H
#define _GCTL_ALGORITHM_H
#include "../lib/core.h"
#include "../lib/io.h"
#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"
using namespace gctl;
#endif // _GCTL_ALGORITHM_H
int main(int argc, char const *argv[]) try
{
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);
}

View File

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

View File

@ -5,12 +5,11 @@ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3")
set(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib)
#
aux_source_directory(geometry/ GCTL_GEOMETRY_SRC)
aux_source_directory(algorithm/ GCTL_ALGORITHM_SRC)
aux_source_directory(core/ GCTL_CORE_SRC)
aux_source_directory(io/ GCTL_IO_SRC)
aux_source_directory(utility/ GCTL_UTILITY_SRC)
aux_source_directory(maths/ GCTL_MATHS_SRC)
aux_source_directory(geodesy/ GCTL_GEODESY_SRC)
aux_source_directory(math/ GCTL_MATH_SRC)
aux_source_directory(graphic/ GCTL_GRAPHIC_SRC)
# windows使getopt_winGNU getopt
# linux
if(NOT WIN32)
@ -24,11 +23,11 @@ endif()
#
#
# libcmake
add_library(gctl SHARED ${GCTL_UTILITY_SRC} ${GCTL_ALGORITHM_SRC} ${GCTL_GEODESY_SRC}
${GCTL_GEOMETRY_SRC} ${GCTL_IO_SRC} ${GCTL_MATHS_SRC})
add_library(gctl SHARED ${GCTL_CORE_SRC} ${GCTL_UTILITY_SRC}
${GCTL_GRAPHIC_SRC} ${GCTL_IO_SRC} ${GCTL_MATH_SRC})
#
add_library(gctl_static STATIC ${GCTL_UTILITY_SRC} ${GCTL_ALGORITHM_SRC} ${GCTL_GEODESY_SRC}
${GCTL_GEOMETRY_SRC} ${GCTL_IO_SRC} ${GCTL_MATHS_SRC})
add_library(gctl_static STATIC ${GCTL_CORE_SRC} ${GCTL_UTILITY_SRC}
${GCTL_GRAPHIC_SRC} ${GCTL_IO_SRC} ${GCTL_MATH_SRC})
#
set_target_properties(gctl_static PROPERTIES OUTPUT_NAME "gctl")
#
@ -42,19 +41,29 @@ 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
find_library(LIB_TAR ${FFTW3_LIB} HINTS ${FFTW3_LIB_DIR})
target_link_libraries(gctl PUBLIC ${LIB_TAR})
target_link_libraries(gctl_static ${LIB_TAR})
find_library(FFTW_LIB ${FFTW3_LIBRARIES} HINTS ${FFTW3_LIBRARY_DIRS})
target_link_libraries(gctl PUBLIC ${FFTW_LIB})
target_link_libraries(gctl_static ${FFTW_LIB})
endif()
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)
@ -102,19 +111,17 @@ endif()
file(GLOB GCTL_HEAD *.h)
file(GLOB GCTL_CORE_HEAD core/*.h)
file(GLOB GCTL_IO_HEAD io/*.h)
file(GLOB GCTL_ALGORITHM_HEAD algorithm/*.h)
file(GLOB GCTL_MATHS_HEAD maths/*.h)
file(GLOB GCTL_GEOMETRY_HEAD geometry/*.h)
file(GLOB GCTL_MATH_HEAD math/*.h)
file(GLOB GCTL_POLY_HEAD poly/*.h)
file(GLOB GCTL_UTILITY_HEAD utility/*.h)
file(GLOB GCTL_MESH_HEAD mesh/*.h)
file(GLOB GCTL_GEODESY_HEAD geodesy/*.h)
file(GLOB GCTL_GRAPHIC_HEAD graphic/*.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_ALGORITHM_HEAD "${PROJECT_SOURCE_DIR}/lib/algorithm/fft.h")
list(REMOVE_ITEM GCTL_MATH_HEAD "${PROJECT_SOURCE_DIR}/lib/math/fft.h")
endif()
if(NOT GCTL_NETCDF)
@ -124,9 +131,7 @@ 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_MATHS_HEAD} DESTINATION include/gctl/maths)
install(FILES ${GCTL_ALGORITHM_HEAD} DESTINATION include/gctl/algorithm)
install(FILES ${GCTL_GEOMETRY_HEAD} DESTINATION include/gctl/geometry)
install(FILES ${GCTL_MATH_HEAD} DESTINATION include/gctl/math)
install(FILES ${GCTL_POLY_HEAD} DESTINATION include/gctl/poly)
install(FILES ${GCTL_UTILITY_HEAD} DESTINATION include/gctl/utility)
install(FILES ${GCTL_MESH_HEAD} DESTINATION include/gctl/mesh)
install(FILES ${GCTL_GEODESY_HEAD} DESTINATION include/gctl/geodesy)
install(FILES ${GCTL_GRAPHIC_HEAD} DESTINATION include/gctl/graphic)

View File

@ -1,139 +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_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 14使
*/
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 14使
*/
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 14使
*/
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

View File

@ -30,12 +30,13 @@
#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

View File

@ -28,14 +28,15 @@
#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*/
@ -62,8 +63,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<double>> _1cd_array; ///< 1-D complex double 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
/**
* @brief Iterator template. This could be used to enable the use of iteration algorithms for array-like objects.
@ -107,7 +108,7 @@ namespace gctl
*
* @param[in] len Length of the array. Must be equal to or bigger than zero.
*/
array(size_t len);
explicit array(size_t len);
/**
* @brief Construct an array with the given length and initial values.
@ -145,7 +146,7 @@ namespace gctl
*
* @param init_val Initial values
*/
array(std::initializer_list<ArrValType> init_val);
explicit array(std::initializer_list<ArrValType> init_val);
/**
* @brief Construct a new array object as a sequence
@ -584,12 +585,6 @@ 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
*
@ -779,7 +774,7 @@ namespace gctl
*
* @return
*/
ArrValType module(norm_type_e n_type = L2);
ArrValType module(norm_type_e n_type = L2) const;
/**
* @brief Normalize the array to the given module value
@ -1509,18 +1504,6 @@ 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)
{
@ -1640,7 +1623,7 @@ namespace gctl
{
if (rn*cn != length_)
{
throw invalid_argument("[gctl::array<T>::sequence2d] Invalid sequence sizes.");
throw std::invalid_argument("[gctl::array<T>::sequence2d] Invalid sequence sizes.");
}
array<ArrValType> out_x(cn), out_y(rn);
@ -1664,7 +1647,7 @@ namespace gctl
{
if (ln*rn*cn != length_)
{
throw invalid_argument("[gctl::array<T>::sequence3d] Invalid sequence sizes.");
throw std::invalid_argument("[gctl::array<T>::sequence3d] Invalid sequence sizes.");
}
array<ArrValType> out_x(cn), out_y(rn), out_z(ln);
@ -1709,7 +1692,7 @@ namespace gctl
throw std::runtime_error("Incompatible array sizes. gctl::array<T>::dot(...)");
}
ArrValType s = 0;
ArrValType s = (ArrValType) 0;
for (size_t i = 0; i < length_; i++)
{
s += val_[i]*a[i];
@ -1744,15 +1727,32 @@ 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 || length_ == 1) return ArrValType{};
if (length_ == 0 || val_ == nullptr) return ArrValType{};
ArrValType m = mean(), s = 0;
ArrValType m = mean(), s = ArrValType(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)
ArrValType array<ArrValType>::module(norm_type_e n_type) const
{
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 runtime_error("[gctl::array<T>::orth] Incompatible array size.");
throw std::runtime_error("[gctl::array<T>::orth] Incompatible array size.");
}
ArrValType product = dot(b);
@ -1913,28 +1913,13 @@ namespace gctl
static_assert(std::is_arithmetic<ArrValType>::value,
"gctl::array<T>::linear2log(...) could only be used with an arithmetic type.");
for (size_t i = 0; i < length_; ++i)
{
val_[i] = log(val_[i])/log(base);
}
return;
}
if (length_ == 0 || val_ == nullptr) 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++)
{
d = d + (val_[i] - mn)*(val_[i] - mn);
if (val_[i] <= ArrValType(0)) continue; // 跳过非正值
val_[i] = std::log(val_[i])/std::log(base);
}
return d/static_cast<ArrValType>(length_);
}
template <typename ArrValType>
@ -1946,7 +1931,7 @@ namespace gctl
if (length_ < 2) return 0;
ArrValType mn = mean();
ArrValType d = 0;
ArrValType d = (ArrValType) 0;
for (size_t i = 0; i < length_; i++)
{
d = d + (val_[i] - mn)*(val_[i] - mn);
@ -2044,6 +2029,6 @@ namespace gctl
else os << byte << " B\n";
return;
}
}
};
#endif // _GCTL_ARRAY_H

View File

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

View File

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

View File

@ -43,12 +43,13 @@
/**
* physical constants
*/
#define GCTL_WGS84_PoleRadius 6356752.314 ///< WGS84椭球极半径
#define GCTL_WGS84_EquatorRadius 6378137 ///< WGS84椭球长半径
#define GCTL_WGS84_FLAT 3.35281066474748e-3 ///< (1/298.25722356300)
#define GCTL_WGS84_PoleRadius 6356752.3142 ///< WGS84椭球极半径
#define GCTL_WGS84_EquatorRadius 6378136.460 ///< WGS84椭球长半径
#define GCTL_WGS84_FLAT 3.35281066474748e-3 ///< (1.0/298.257223563)
#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

View File

@ -28,25 +28,22 @@
#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 array types.
* Here we define some commonly used matrix types.
*/
typedef matrix<int> _2i_matrix; ///< 2-D integer array
typedef matrix<float> _2f_matrix; ///< 2-D single-precision floating-point array
typedef matrix<double> _2d_matrix; ///< 2-D double-precision floating-point array
typedef matrix<std::string> _2s_matrix; ///< 2-D string array
typedef matrix<std::complex<float>> _2cf_matrix; ///< 1-D complex float array
typedef matrix<std::complex<double>> _2cd_matrix; ///< 1-D complex double array
typedef matrix<int> _2i_matrix; ///< 2-D integer matrix
typedef matrix<float> _2f_matrix; ///< 2-D single-precision floating-point matrix
typedef matrix<double> _2d_matrix; ///< 2-D double-precision floating-point matrix
typedef matrix<std::string> _2s_matrix; ///< 2-D string matrix
typedef matrix<std::complex<float> > _2cf_matrix; ///< 1-D complex float matrix
typedef matrix<std::complex<double> > _2cd_matrix; ///< 1-D complex double matrix
/**
* @brief
@ -117,7 +114,7 @@ namespace gctl
*
* @param[in] b
*/
matrix(const std::vector<std::vector<MatValType> > &b);
explicit matrix(const std::vector<std::vector<MatValType> > &b);
/**
* @brief
@ -248,6 +245,22 @@ 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
*
@ -263,12 +276,19 @@ namespace gctl
*/
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 norm(norm_type_e n_type = L2) const;
double module(norm_type_e n_type = L2) const;
/**
* @brief Set elements' value as a sequent.
@ -786,6 +806,32 @@ 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)
{
@ -832,7 +878,47 @@ namespace gctl
}
template <typename MatValType>
double matrix<MatValType>::norm(norm_type_e n_type) const
void matrix<MatValType>::normalize(matrix_order_e odr)
{
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)
@ -1096,6 +1182,6 @@ namespace gctl
}
return;
}
}
};
#endif // _GCTL_MATRIX_H

View File

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

View File

@ -28,8 +28,6 @@
#ifndef _GCTL_SPMAT_H
#define _GCTL_SPMAT_H
#include "enum.h"
#include "array.h"
#include "matrix.h"
namespace gctl
@ -42,8 +40,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<double>> _1cd_spmat; ///< 双精度复浮点稀疏矩阵
typedef spmat<std::complex<float> > _1cf_spmat; ///< 单精度复浮点稀疏矩阵
typedef spmat<std::complex<double> > _1cd_spmat; ///< 双精度复浮点稀疏矩阵
/**
* @brief
@ -63,6 +61,7 @@ namespace gctl
*/
mat_node()
{
val = (NodeValueType) 0;
r_id = c_id = -1;
next_right = next_down = nullptr;
}
@ -628,6 +627,7 @@ namespace gctl
r_end = c_end = nullptr;
d_link = nullptr;
r_count = c_count = nullptr;
zero_val = (MatrixValueType) 0;
}
template <typename MatrixValueType>
@ -2941,6 +2941,7 @@ 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];
@ -3304,6 +3305,6 @@ namespace gctl
#endif // GCTL_EIGEN
}
};
#endif //_GCTL_SPMAT_H

View File

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

207
lib/core/str.cpp Normal file
View File

@ -0,0 +1,207 @@
/********************************************************
*
*
*
*
*
*
* 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;
}

View File

@ -1,4 +1,4 @@
/********************************************************
/********************************************************
*
*
*
@ -25,32 +25,18 @@
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#ifndef _GCTL_STREAM_TEMPLATE_H
#define _GCTL_STREAM_TEMPLATE_H
#ifndef _GCTL_STR_H
#define _GCTL_STR_H
#include "string"
#include "sstream"
#include "exceptions.h"
#include "../core/exceptions.h"
#include "../core/vector_t.h"
#include <cmath>
#include <cstring>
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
/**
* terminal control symbols
*/
#define GCTL_BOLD "\033[1m" ///< 设置终端字体与颜色为粗体
#define GCTL_BOLDRED "\033[1m\033[31m" ///< 设置终端字体与颜色为粗体红色
#define GCTL_BOLDGREEN "\033[1m\033[32m" ///< 设置终端字体与颜色为粗体绿色
#define GCTL_BOLDBLUE "\033[1m\033[34m" ///< 设置终端字体与颜色为粗体蓝色
#define GCTL_BOLDYELLOW "\033[1m\033[33m" ///< 设置后续字符字体为黄色加粗
#define GCTL_UNDERLINE "\033[1m\033[4m" ///< 设置终端字体为下划线
#define GCTL_RESET "\033[0m" ///< 重置字体与颜色设定
#define GCTL_CLEARLINE "\033[K" ///< 清空终端当前行
#define GCTL_CLEARALL "\033[2J" ///< 清空终端窗口
#define GCTL_MOVEUP(os, x) do {os << "\033[" << x << "A";} while (0); ///< 终端光标上移x行
#define GCTL_MOVEDOWN(os, x) do {os << "\033[" << x << "B";} while (0); ///< 终端光标下移x行
#define GCTL_MOVERIGHT(os, x) do {os << "\033[" << x << "C";} while (0); ///< 终端光标右移x列
#define GCTL_MOVELEFT(os, x) do {os << "\033[" << x << "D";} while (0); ///< 终端光标左移x列
#define GCTL_MOVETO(os, r, c) do {os << "\033[" << c << ";" << r << "H";} while (0); ///< 终端光标移动至r行c列位置
namespace gctl
{
template <typename ObjValType>
@ -132,7 +118,7 @@ namespace gctl
}
template <typename ValueType>
void str2type_vector2d(const std::vector<std::vector<std::string>> &vec2d,
void str2type_vector2d(const std::vector<std::vector<std::string> > &vec2d,
std::string order_str, std::vector<ValueType> &vec)
{
int idx;
@ -174,7 +160,7 @@ namespace gctl
}
template <typename ValueType>
void type2str_vector2d(std::vector<std::vector<ValueType>> &vec2d,
void type2str_vector2d(std::vector<std::vector<ValueType> > &vec2d,
const std::vector<std::string> &vec)
{
std::vector<std::string> str_vec;
@ -185,7 +171,7 @@ namespace gctl
}
template <typename ValueType, typename... Args>
void type2str_vector2d(std::vector<std::vector<ValueType>> &vec2d,
void type2str_vector2d(std::vector<std::vector<ValueType> > &vec2d,
const std::vector<std::string> &vec, Args&... rest)
{
std::vector<std::string> str_vec;
@ -197,52 +183,6 @@ namespace gctl
return;
}
/**
* @brief Parse a string argument into an element array
*
* @param[in] val_str The input value string
* @param out_vec The output array
* @param[in] separator The separator
*
* @tparam ValueType Runtime value type
*
* @return 0 for success, -1 for failure
*/
template <typename ValueType>
int parse_string_to_vector(std::string val_str, char separator, std::vector<ValueType> &out_vec)
{
val_str.erase(0, val_str.find_first_not_of(" \t"));
val_str.erase(val_str.find_last_not_of(" \t") + 1);
int pos;
bool bk = false;
ValueType val;
std::string l_str;
while (1)
{
pos = val_str.find(separator);
if (pos == val_str.npos)
{
l_str = val_str;
bk = true;
}
else
{
l_str = val_str.substr(0, pos);
val_str = val_str.substr(pos+1, val_str.length());
val_str.erase(0, val_str.find_first_not_of(" \t"));
}
// 如果输出也是字符串 就直接赋值即可避免l_str中含有空格会出现bug
if constexpr (std::is_same<ValueType, std::string>::value) val = l_str;
else str2type(l_str, val);
out_vec.push_back(val);
if (bk) break;
}
return 0;
}
/**
* @brief Parse a string argument into separated elements.
*
@ -317,6 +257,132 @@ namespace gctl
return parse_string_to_value(tmp_str, separator, throw_err, rest...) + 1;
}
}
#endif // _GCTL_STREAM_TEMPLATE_H
/**
* @brief
*
* @param new_str
* @param old_str
* @param old_value
* @param new_value
*
* @return
*/
int replace_all(std::string& new_str, const std::string& old_str, const std::string& old_value,
const std::string& new_value);
/**
* @brief
*
* @param in_str
* @param patch_str
* @return
*/
std::string patch_string(std::string in_str, std::string patch_str);
/**
* @brief string对象为stringstream对象
*
* @param[in] in_str string字符串
* @param out_ss stringstreams对象
* @param[in] delimiter string对象时的分隔符
* delimiter转换为空格
*/
void str2ss(std::string in_str, std::stringstream &out_ss, std::string delimiter = "");
/**
* @brief string字符串为double类型的数值
*
* nan或者inf等表示无效值的符号
* >>
* fortran输出文件中以D标注的科学数字转换为浮点类型
*
* @param[in] instr
*
* @return double类型的数值
*/
double str2double(std::string instr);
/**
* @brief double类型数值为string类型字符串 to_string函数
*
* @param[in] in_d
*
* @return
*/
std::string double2str(double in_d);
/**
* @brief
*
* @param[in] length
* @param out
*/
void random_char(unsigned int length, char* out);
/**
* @brief
*
* @param[in] length
* @param out_str
*/
void random_str(unsigned int lenght, std::string &out_str);
/**
* @brief Parse a string argument into an element array
*
* @param[in] val_str The input value string
* @param out_vec The output array
* @param[in] separator The separator
*
* @tparam ValueType Runtime value type
*
* @return 0 for success, -1 for failure
*/
template <typename ValueType>
int parse_string_to_vector(std::string val_str, char separator, std::vector<ValueType> &out_vec)
{
val_str.erase(0, val_str.find_first_not_of(" \t"));
val_str.erase(val_str.find_last_not_of(" \t") + 1);
int pos;
bool bk = false;
ValueType val;
std::string l_str;
while (1)
{
pos = val_str.find(separator);
if (pos == val_str.npos)
{
l_str = val_str;
bk = true;
}
else
{
l_str = val_str.substr(0, pos);
val_str = val_str.substr(pos+1, val_str.length());
val_str.erase(0, val_str.find_first_not_of(" \t"));
}
// 如果输出也是字符串 就直接赋值即可避免l_str中含有空格会出现bug
if constexpr (std::is_same<ValueType, std::string>::value) val = l_str;
else str2type(l_str, val);
out_vec.push_back(val);
if (bk) break;
}
return 0;
}
/**
* @brief
*
* @note
*
* @param in_str
* @param str_vec
*/
void parse_string_with_quotes(std::string in_str, std::vector<std::string> &str_vec);
};
#endif // _GCTL_STR_H

View File

@ -28,9 +28,9 @@
#ifndef _GCTL_VECTOR_TYPE_H
#define _GCTL_VECTOR_TYPE_H
#include "vector"
#include "complex"
#include "iostream"
#include <vector>
#include <complex>
#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<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<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<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<double>>> _2cd_vector; ///< 2-D complex double 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::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<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::complex<float> > > _2cf_vector; ///< 2-D complex float 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

View File

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

283
lib/graphic/cliplot.cpp Normal file
View File

@ -0,0 +1,283 @@
/********************************************************
*
*
*
*
*
*
* 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;
}

206
lib/graphic/cliplot.h Normal file
View File

@ -0,0 +1,206 @@
/********************************************************
*
*
*
*
*
*
* 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

148
lib/graphic/gmt.cpp Normal file
View File

@ -0,0 +1,148 @@
/********************************************************
*
*
*
*
*
*
* 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;
}

119
lib/graphic/gmt.h Normal file
View File

@ -0,0 +1,119 @@
/********************************************************
*
*
*
*
*
*
* 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

135
lib/graphic/gnuplot.cpp Normal file
View File

@ -0,0 +1,135 @@
/********************************************************
*
*
*
*
*
*
* 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;
}

103
lib/graphic/gnuplot.h Normal file
View File

@ -0,0 +1,103 @@
/********************************************************
*
*
*
*
*
*
* 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

View File

@ -28,15 +28,21 @@
#ifndef _GCTL_IO_H
#define _GCTL_IO_H
#include "io/mesh_io.h"
#include "io/term_io.h"
#include "io/file_io.h"
#include "io/native_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"
#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"
#endif // _GCTL_IO_H

View File

@ -29,7 +29,6 @@
gctl::dsv_io::dsv_io()
{
file_ = "";
att_sym_ = '#';
tag_sym_ = '!';
deli_sym_ = ' ';
@ -43,9 +42,8 @@ gctl::dsv_io::~dsv_io()
clear();
}
gctl::dsv_io::dsv_io(std::string filename, std::string file_exten, table_headtype_e t)
gctl::dsv_io::dsv_io(std::string filename, std::string file_exten, int t)
{
file_ = "";
att_sym_ = '#';
tag_sym_ = '!';
deli_sym_ = ' ';
@ -59,7 +57,6 @@ gctl::dsv_io::dsv_io(std::string filename, std::string file_exten, table_headtyp
void gctl::dsv_io::clear()
{
file_ = "";
att_sym_ = '#';
tag_sym_ = '!';
deli_sym_ = ' ';
@ -74,27 +71,17 @@ void gctl::dsv_io::clear()
return;
}
void gctl::dsv_io::get_row_names(std::vector<std::string> &names)
std::vector<std::string> gctl::dsv_io::row_names()
{
names.resize(row_num_);
_1s_vector names(row_num_);
for (size_t i = 1; i <= row_num_; i++)
{
names[i - 1] = table_[i][0].str_;
}
return;
return names;
}
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)
void gctl::dsv_io::row_names(const std::vector<std::string> &names, const std::vector<int> &idx, std::string corner_name)
{
if (!idx.empty())
{
@ -117,7 +104,17 @@ void gctl::dsv_io::set_row_names(const std::vector<std::string> &names, const st
return;
}
void gctl::dsv_io::set_column_names(const std::vector<std::string> &names, const std::vector<int> &idx)
std::vector<std::string> gctl::dsv_io::column_names()
{
_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())
{
@ -138,11 +135,11 @@ void gctl::dsv_io::set_column_names(const std::vector<std::string> &names, const
return;
}
void gctl::dsv_io::set_row_type(cell_type_e t, int idx)
void gctl::dsv_io::row_type(table_cell_type t, int idx)
{
if (idx > row_num_ || idx <= 0)
{
throw std::runtime_error("[gctl::dsv_io] Invalid column index.");
throw std::runtime_error("[gctl::dsv_io] Invalid row index.");
}
for (size_t i = 0; i <= col_num_; i++)
@ -151,12 +148,27 @@ void gctl::dsv_io::set_row_type(cell_type_e t, int idx)
}
}
void gctl::dsv_io::set_row_type(cell_type_e t, std::string name)
void gctl::dsv_io::row_type(table_cell_type t, std::string name)
{
set_row_type(t, name_index(name, true));
row_type(t, name_index(name, true));
}
void gctl::dsv_io::set_column_type(cell_type_e t, int idx)
gctl::table_cell_type gctl::dsv_io::row_type(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)
{
@ -170,39 +182,53 @@ void gctl::dsv_io::set_column_type(cell_type_e t, int idx)
return;
}
void gctl::dsv_io::set_column_type(cell_type_e t, std::string name)
void gctl::dsv_io::column_type(table_cell_type t, std::string name)
{
set_column_type(t, name_index(name, false));
column_type(t, name_index(name, false));
}
void gctl::dsv_io::load_text(std::string filename, std::string file_exten, table_headtype_e t)
gctl::table_cell_type gctl::dsv_io::column_type(int idx)
{
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_) //读入到头信息中
@ -225,7 +251,7 @@ void gctl::dsv_io::load_text(std::string filename, std::string file_exten, table
// 动态调整列数
cn = tmp_cols.size();
cn_max = std::max(cn, cn_max);
cn_max = GCTL_MAX(cn, cn_max);
table_[i].resize(tmp_cols.size());
for (size_t j = 0; j < tmp_cols.size(); j++)
@ -260,7 +286,7 @@ void gctl::dsv_io::load_text(std::string filename, std::string file_exten, table
}
}
if (t == ColumnHead) // 有列头 需要补齐空白的行头
if (t == ColHead) // 有列头 需要补齐空白的行头
{
row_num_ = table_.size() - 1;
col_num_ = table_[0].size();
@ -280,7 +306,7 @@ void gctl::dsv_io::load_text(std::string filename, std::string file_exten, table
table_[0].resize(col_num_ + 1);
}
if (t == BothHead) // 有行头和列头
if ((t & RowHead) && (t & ColHead)) // 有行头和列头
{
row_num_ = table_.size() - 1;
col_num_ = table_[0].size() - 1;
@ -290,9 +316,9 @@ void gctl::dsv_io::load_text(std::string filename, std::string file_exten, table
return;
}
void gctl::dsv_io::load_csv(std::string filename, table_headtype_e t)
void gctl::dsv_io::load_csv(std::string filename, int t)
{
set_delimeter(',');
delimeter(',');
load_text(filename, ".csv", t);
return;
}
@ -317,80 +343,30 @@ 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++)
{
if (table_[i][0].out_ok_ && table_[i][0].str_ != "")
for (size_t j = 0; j <= col_num_; j++)
{
col_st = 0;
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
if (table_[i][j].out_ok_ && table_[i][j].str_!= "")
{
outfile << table_[i][j].str_;
line_st = true; // start line
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;
}
}
outfile << std::endl;
}
*/
outfile.close();
return;
}
void gctl::dsv_io::save_csv(std::string filename)
{
set_delimeter(',');
delimeter(',');
save_text(filename, ".csv");
return;
}
@ -411,72 +387,167 @@ void gctl::dsv_io::init_table(int row, int col)
return;
}
void gctl::dsv_io::info(table_headtype_e t)
gctl::dsv_io gctl::dsv_io::export_table(bool ignore_disabled)
{
std::clog << "File: " << file_ << "\n------------\n";
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::vector<std::string> str_line, row_names, col_names;
std::vector<std::vector<std::string> > str_table;
if (t == ColumnHead || t == BothHead)
std::string cor_name = table_[0][0].str_;
for (size_t j = 1; j <= col_num_; j++)
{
std::clog << "Columns:\n";
if (table_[0][j].out_ok_ || !ignore_disabled)
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_ != "")
{
std::clog << table_[0][i].str_ << " | ";
if (table_[1][i].type_ == String) std::clog << "String | ";
if (table_[1][i].type_ == Int) std::clog << "Int | ";
if (table_[1][i].type_ == Float) std::clog << "Float | ";
std::clog << table_[1][i].str_ << " -> " << table_[row_num_][i].str_;
os << table_[0][i].str_ << " | ";
if (table_[0][i].out_ok_) os << "Enabled | ";
else os << "Disabled | ";
if (table_[1][i].type_ == String) os << "String | ";
if (table_[1][i].type_ == Int) os << "Int | ";
if (table_[1][i].type_ == Float) os << "Float | ";
os << table_[1][i].str_ << " -> " << table_[row_num_][i].str_;
}
else
{
std::clog << "C" + std::to_string(i) << " | ";
if (table_[1][i].type_ == String) std::clog << "String | ";
if (table_[1][i].type_ == Int) std::clog << "Int | ";
if (table_[1][i].type_ == Float) std::clog << "Float | ";
std::clog << table_[1][i].str_ << " -> " << table_[row_num_][i].str_;
os << "C" + std::to_string(i) << " | ";
if (table_[0][i].out_ok_) os << "Enabled | ";
else os << "Disabled | ";
if (table_[1][i].type_ == String) os << "String | ";
if (table_[1][i].type_ == Int) os << "Int | ";
if (table_[1][i].type_ == Float) os << "Float | ";
os << table_[1][i].str_ << " -> " << table_[row_num_][i].str_;
}
if (!table_[0][i].out_ok_) std::clog << " (No output)";
std::clog << std::endl;
os << std::endl;
}
std::clog << "============\n";
os << "------------\n";
}
if (t == RowHead || t == BothHead)
if (t & RowInfo)
{
std::clog << "Rows:\n";
os << "Rows:\n";
for (size_t i = 1; i <= row_num_; i++)
{
if (table_[i][0].str_ != "")
{
std::clog << table_[i][0].str_ << " | ";
if (table_[i][1].type_ == String) std::clog << "String | ";
if (table_[i][1].type_ == Int) std::clog << "Int | ";
if (table_[i][1].type_ == Float) std::clog << "Float | ";
std::clog << table_[i][1].str_ << " -> " << table_[i][col_num_].str_;
os << table_[i][0].str_ << " | ";
if (table_[i][0].out_ok_) os << "Enabled | ";
else os << "Disabled | ";
if (table_[i][1].type_ == String) os << "String | ";
if (table_[i][1].type_ == Int) os << "Int | ";
if (table_[i][1].type_ == Float) os << "Float | ";
os << table_[i][1].str_ << " -> " << table_[i][col_num_].str_;
}
else
{
std::clog << "R" + std::to_string(i) << " | ";
if (table_[i][1].type_ == String) std::clog << "String | ";
if (table_[i][1].type_ == Int) std::clog << "Int | ";
if (table_[i][1].type_ == Float) std::clog << "Float | ";
std::clog << table_[i][1].str_ << " -> " << table_[i][col_num_].str_;
os << "R" + std::to_string(i) << " | ";
if (table_[i][0].out_ok_) os << "Enabled | ";
else os << "Disabled | ";
if (table_[i][1].type_ == String) os << "String | ";
if (table_[i][1].type_ == Int) os << "Int | ";
if (table_[i][1].type_ == Float) os << "Float | ";
os << table_[i][1].str_ << " -> " << table_[i][col_num_].str_;
}
if (!table_[i][0].out_ok_) std::clog << " (No output)";
std::clog << std::endl;
os << std::endl;
}
std::clog << "============\n";
os << "------------\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>
@ -517,6 +588,19 @@ 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)
@ -623,162 +707,63 @@ 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::filt_column(std::string cnd_str, std::string cnd_col,
const std::vector<std::string> &out_col, dsv_io &out_table)
void gctl::dsv_io::filter(std::string cnd_str, std::string cnd_tar, table_headtype_e thead)
{
int idx = name_index(cnd_col);
if (idx < 0) throw std::runtime_error("[gctl::dsv_io::] Invalid column index or name.");
int idx;
if (thead == RowHead) idx = name_index(cnd_tar, true);
else if (thead == ColHead) idx = name_index(cnd_tar);
else throw std::runtime_error("[gctl::dsv_io::filter] Invalid table head type.");
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.");
}
}
if (idx < 0) throw std::runtime_error("[gctl::dsv_io::filter] Invalid row/column index or name.");
std::smatch ret;
std::regex pat(cnd_str);
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++)
if (thead == RowHead) // cnd_tar是行头 此时为按列过滤
{
if (regex_search(table_[i][idx].str_, ret, pat))
for (size_t i = 1; i <= col_num_; i++)
{
if (out_row)
if (!regex_search(table_[idx][i].str_, ret, pat))
{
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);
column_output(i, Disable);
}
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_);
}
}
out_table.init_table(str_table);
std::vector<std::string> io_col;
if (out_row)
else // cnd_tar是列头 此时为按行过滤
{
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++)
for (size_t i = 1; i <= row_num_; i++)
{
io_col.push_back(out_col[j]);
if (!regex_search(table_[i][idx].str_, ret, pat))
{
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::filt_column(rowbool_func_t func, const std::vector<std::string> &out_col, dsv_io &out_table)
void gctl::dsv_io::filter(linebool_func_t func, table_headtype_e thead)
{
array<int> odx;
bool out_row = false;
if (out_col.empty()) out_row = true;
else
if (thead == RowHead)
{
odx.resize(out_col.size());
for (size_t i = 0; i < out_col.size(); i++)
for (size_t i = 1; i <= row_num_; 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 (!func(table_[i])) row_output(i, Disable);
}
}
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++)
else if (thead == ColHead)
{
if (func(table_[i]))
std::vector<table_cell> col_cell(row_num_);
for (size_t i = 1; i <= col_num_; i++)
{
if (out_row)
for (size_t j = 1; j < row_num_; j++)
{
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 < odx.size(); j++)
{
str_line.push_back(table_[i][odx[j]].str_);
}
str_table.push_back(str_line);
col_cell[j] = table_[j][i];
}
row_names.push_back(table_[i][0].str_);
if (!func(col_cell)) column_output(i, Disable);
}
}
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);
else throw std::runtime_error("[gctl::dsv_io::filter] Invalid table head type.");
return;
}
@ -830,41 +815,45 @@ void gctl::dsv_io::cal_column(std::string expr_str, const std::vector<std::strin
return;
}
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)
void gctl::dsv_io::filter(std::string cnd_str, const std::vector<std::string> &cnd_tars, table_headtype_e thead)
{
array<int> idx(cnd_col.size());
for (size_t i = 0; i < cnd_col.size(); i++)
array<int> idx(cnd_tars.size());
if (thead == RowHead)
{
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)
for (size_t i = 0; i < cnd_tars.size(); i++)
{
throw std::runtime_error("[gctl::dsv_io] Invalid column type for numerical calculating.");
idx[i] = name_index(cnd_tars[i], true);
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.");
}
}
}
array<int> odx;
bool out_row = false;
if (out_col.empty()) out_row = true;
else
else if (thead == ColHead)
{
odx.resize(out_col.size());
for (size_t i = 0; i < out_col.size(); i++)
for (size_t i = 0; i < cnd_tars.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.");
idx[i] = name_index(cnd_tars[i]);
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_col.size());
array<double> var(cnd_tars.size());
for (size_t i = 0; i < var.size(); i++)
{
symbol_table.add_variable(cnd_col[i], var[i]);
symbol_table.add_variable(cnd_tars[i], var[i]);
}
exprtk::expression<double> expression;
@ -876,77 +865,32 @@ void gctl::dsv_io::filt_column(std::string cnd_str, const std::vector<std::strin
throw std::runtime_error("[gctl::dsv_io] Fail to compile the math expression.");
}
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++)
if (thead == RowHead) // cnd_tars是行头 此时为按列过滤
{
for (size_t j = 0; j < var.size(); j++)
for (size_t i = 1; i <= col_num_; i++)
{
var[j] = table_[i][idx[j]].value<double>();
}
if (expression.value() > 0.5) // return 1 if matched or 0 if dismatched
{
if (out_row)
for (size_t j = 0; j < var.size(); j++)
{
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);
var[j] = table_[idx[j]][i].value<double>();
}
row_names.push_back(table_[i][0].str_);
// return 1 if matched or 0 if dismatched
if (expression.value() < 0.5) column_output(i, Disable);
}
}
out_table.init_table(str_table);
std::vector<std::string> io_col;
if (out_row)
else // cnd_tars是列头 此时为按行过滤
{
get_column_names(io_col);
out_table.cell(table_[0][0].str_, 0, 0);
}
else
{
for (size_t j = 0; j < idx.size(); j++)
for (size_t i = 1; i <= row_num_; i++)
{
io_col.push_back(cnd_col[j]);
}
for (size_t j = 0; j < var.size(); j++)
{
var[j] = table_[i][idx[j]].value<double>();
}
for (size_t j = 0; j < odx.size(); j++)
{
io_col.push_back(out_col[j]);
// return 1 if matched or 0 if dismatched
if (expression.value() < 0.5) 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;
}
@ -956,9 +900,8 @@ gctl::geodsv_io::geodsv_io(){}
gctl::geodsv_io::~geodsv_io(){}
gctl::geodsv_io::geodsv_io(std::string filename, std::string file_exten, table_headtype_e t)
gctl::geodsv_io::geodsv_io(std::string filename, std::string file_exten, int t)
{
file_ = "";
att_sym_ = '#';
tag_sym_ = '!';
deli_sym_ = ' ';

View File

@ -28,19 +28,23 @@
#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 cell_type_e
enum table_cell_type
{
String,
Int,
@ -50,7 +54,7 @@ namespace gctl
struct table_cell
{
std::string str_; // 单元格的内容 统一保存为字符串
cell_type_e type_; // 类型字符串
table_cell_type type_; // 类型字符串
bool out_ok_; // 是否可输出到文件
table_cell()
@ -101,15 +105,25 @@ namespace gctl
};
/**
* @brief
*
* @brief
*/
enum table_headtype_e
{
NoHead, // 没有表头
BothHead, // 同时有行与列表头
ColumnHead, // 只有列表头
RowHead, // 只有行表头
NoHead = 1, // 0001 无表头
ColHead = 2, // 0010 有列表头
RowHead = 4, // 0100 有行表头
};
/**
* @brief
*/
enum table_infotype_e
{
TagInfo = 1, // 00001
AttInfo = 2, // 00010
HeadInfo = 4, // 00100
RowInfo = 8, // 01000
ColInfo = 16, // 10000
};
/**
@ -126,8 +140,6 @@ namespace gctl
class dsv_io
{
protected:
// 文件名
std::string file_;
// 头信息行数 表格行数(不包括表头) 表格列数(不包括表头)
int head_num_, row_num_, col_num_;
// 注释行起始符 标记行起始符 分割符
@ -157,21 +169,13 @@ namespace gctl
* @param file_exten
* @param t
*/
dsv_io(std::string filename, std::string file_exten = ".txt", table_headtype_e t = NoHead);
dsv_io(std::string filename, std::string file_exten = ".txt", int t = NoHead);
/**
* @brief
*
* @brief
*/
void clear();
/**
* @brief
*
* @return
*/
int head_number(){return head_num_;}
/**
* @brief
*
@ -187,88 +191,102 @@ namespace gctl
int col_number(){return col_num_;}
/**
* @brief
* @brief
*
* @return
* @return
*/
const std::vector<std::string>& get_head_records(){return heads_;}
/**
* @brief
*
* @return
*/
const std::vector<std::string>& get_annotoations(){return annotates_;}
/**
* @brief
*
* @return
*/
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;}
int head_number(){return head_num_;}
/**
* @brief
*
* @param num
*/
void set_head_number(char num){head_num_ = num;}
void head_number(int num){head_num_ = num;}
/**
* @brief
* @brief
*
* @param att_sym
* @return
*/
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;}
const std::vector<std::string> &head_records(){return heads_;}
/**
* @brief
*
* @param heads
*/
void set_head_records(const std::vector<std::string> &heads){heads_ = heads; head_num_ = heads_.size();}
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 set_annotoations(const std::vector<std::string> &att){annotates_ = att;}
void annotations(const std::vector<std::string> &att){annotates_ = att;}
/**
* @brief
*
* @return
*/
const std::vector<std::string> &annotations(){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 set_tags(const std::vector<std::string> &tags){tags_ = tags;}
void tags(const std::vector<std::string> &tags){tags_ = tags;}
/**
* @brief
*
* @return
*/
const std::vector<std::string> &tags(){return tags_;}
/**
* @brief
@ -277,7 +295,14 @@ namespace gctl
* @param idx
* @param corner_name RowNames
*/
void set_row_names(const std::vector<std::string> &names, const std::vector<int> &idx = {}, std::string corner_name = "RowNames");
void 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
@ -285,7 +310,14 @@ namespace gctl
* @param names
* @param idx
*/
void set_column_names(const std::vector<std::string> &names, const std::vector<int> &idx = {});
void column_names(const std::vector<std::string> &names, const std::vector<int> &idx = {});
/**
* @brief
*
* @return
*/
std::vector<std::string> column_names();
/**
* @brief
@ -293,7 +325,7 @@ namespace gctl
* @param t String|Float|Int
* @param idx
*/
void set_row_type(cell_type_e t, int idx);
void row_type(table_cell_type t, int idx);
/**
* @brief
@ -301,7 +333,23 @@ namespace gctl
* @param t String|Float|Int
* @param name
*/
void set_row_type(cell_type_e t, std::string name);
void row_type(table_cell_type 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
@ -309,7 +357,7 @@ namespace gctl
* @param t String|Float|Int
* @param idx
*/
void set_column_type(cell_type_e t, int idx);
void column_type(table_cell_type t, int idx);
/**
* @brief
@ -317,7 +365,23 @@ namespace gctl
* @param t String|Float|Int
* @param name
*/
void set_column_type(cell_type_e t, std::string name);
void column_type(table_cell_type 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
@ -326,14 +390,14 @@ namespace gctl
* @param file_exten
* @param t
*/
void load_text(std::string filename, std::string file_exten = ".txt", table_headtype_e t = NoHead);
void load_text(std::string filename, std::string file_exten = ".txt", int t = NoHead);
/**
* @brief CSV文件
*
* @param filename
*/
void load_csv(std::string filename, table_headtype_e t = ColumnHead);
void load_csv(std::string filename, int t = ColHead);
/**
* @brief
@ -358,22 +422,43 @@ namespace gctl
*/
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(table_headtype_e t = ColumnHead);
void info(int t, std::ostream &os = std::cout);
/**
* @brief
*
*/
void display();
/**
* @brief name和R<id>C<id>
*
* @param name R<id>C<id>
* @param iter_row name参数为R<id>C<id>
* @param iter_row
* @return 1 -1
*/
int name_index(std::string name, bool iter_row = false);
/**
* @brief
*
* @param s
*/
void table_output(switch_type_e s);
/**
* @brief 使
*
@ -447,32 +532,27 @@ namespace gctl
int add_row(std::string name, std::string id_name);
/**
* @brief
*
* @note
* @brief
*
* @param cnd_str
* @param cnd_col
* @param out_col
* @param out_table
* @param cnd_tar
* @param thead RowHead表示名称为行头则按列过滤ColHead时表示名称为列头则按行过滤
*/
void filt_column(std::string cnd_str, std::string cnd_col,
const std::vector<std::string> &out_col, dsv_io &out_table);
void filter(std::string cnd_str, std::string cnd_tar, table_headtype_e thead = RowHead);
/**
* @brief row operate function pointer
* @brief table line operate function pointer
*
*/
typedef bool (*rowbool_func_t)(const std::vector<table_cell> &table_row);
typedef bool (*linebool_func_t)(const std::vector<table_cell> &table_line);
/**
* @brief
* @brief
*
* @param func
* @param out_col
* @param out_table
* @param func
* @param thead RowHead时表示按行过滤ColHead时表示按列过滤
*/
void filt_column(rowbool_func_t func, const std::vector<std::string> &out_col, dsv_io &out_table);
void filter(linebool_func_t func, table_headtype_e thead = RowHead);
#ifdef GCTL_EXPRTK
@ -489,18 +569,16 @@ 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_col
* @param out_col
* @param out_table
* @param cnd_tars
* @param thead RowHead表示名称为行头则按列过滤ColHead时表示名称为列头则按行过滤
*/
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);
void filter(std::string cnd_str, const std::vector<std::string> &cnd_tars, table_headtype_e thead = RowHead);
#endif // GCTL_EXPRTK
@ -805,7 +883,7 @@ namespace gctl
* @param filename
* @param file_exten
*/
geodsv_io(std::string filename, std::string file_exten = ".txt", table_headtype_e t = NoHead);
geodsv_io(std::string filename, std::string file_exten = ".txt", int t = NoHead);
/**
* @brief
@ -967,6 +1045,6 @@ namespace gctl
*/
void get_column_point3ds(array<point3ds> &data, std::string rname, std::string pname, std::string tname);
};
}
};
#endif //_GCTL_DSV_IO_H

View File

@ -1,4 +1,4 @@
/********************************************************
/********************************************************
*
*
*
@ -25,166 +25,7 @@
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#include "stream.h"
#include "cmath"
#if defined _WINDOWS || __WIN32__
#include "io.h"
// Test for file existence
#define F_OK 0
#endif
//在终端显示一个简易的GCTL图标
void gctl::display_logo(std::ostream &sout)
{
sout << " ___ ___ _____ __\n";
sout << " / _ \\ / __\\/__ \\ / /\n";
sout << " / /_\\// / / /\\// /\n";
sout << "/ /_\\\\/ /___ / / / /___\n";
sout << "\\____/\\____/ \\/ \\____/\n";
sout << "Geophysical Computational Tools & Library\n";
sout << "Author: Dr. Yi Zhang (yizhang-geo@zju.edu.cn)\n\n";
return;
}
//替换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;
}
#include "file_io.h"
//测试打开输入文件 如果成功则返回0并输出信息 否则返回1
void gctl::open_infile(std::ifstream &infile, std::string filename, std::string extension,
@ -192,7 +33,7 @@ void gctl::open_infile(std::ifstream &infile, std::string filename, std::string
{
if (filename == "")
{
throw domain_error("Empty file name. From open_infile(...)");
throw std::domain_error("Empty file name. From open_infile(...)");
}
// 文件流开启则先关闭文件
@ -209,7 +50,7 @@ void gctl::open_infile(std::ifstream &infile, std::string filename, std::string
infile.open(filename+extension, mode);
if (!infile)
{
throw domain_error("Fail to open file: " + filename + extension + ". From open_infile(...)");
throw gctl::domain_error("Fail to open file: " + filename + extension + ". From open_infile(...)");
}
return;
}
@ -219,7 +60,7 @@ void gctl::open_infile(std::ifstream &infile, std::string filename, std::string
infile.open(filename, mode);
if (!infile)
{
throw domain_error("Fail to open file: " + filename + ". From open_infile(...)");
throw gctl::domain_error("Fail to open file: " + filename + ". From open_infile(...)");
}
return;
}
@ -229,7 +70,7 @@ void gctl::open_matched_infile(std::ifstream &infile, std::string filename,
{
if (exten_pattern == "" || filename == "")
{
throw domain_error("Empty file name or extension(s). From open_matched_infile(...)");
throw std::domain_error("Empty file name or extension(s). From open_matched_infile(...)");
}
if (infile.is_open()) infile.close();
@ -354,81 +195,6 @@ void gctl::open_matched_outfile(std::ofstream &outfile, std::string filename,
return;
}
int gctl::terminal_width()
{
int width;
//获取终端窗口的行列数
#ifdef _WINDOWS
CONSOLE_SCREEN_BUFFER_INFO csbi;
GetConsoleScreenBufferInfo(GetStdHandle(STD_OUTPUT_HANDLE), &csbi);
width = csbi.srWindow.Right - csbi.srWindow.Left + 1;
#else
struct winsize w;
ioctl(STDOUT_FILENO, TIOCGWINSZ, &w);
width = w.ws_col;
#endif
return width;
}
int gctl::terminal_height()
{
int height;
//获取终端窗口的行列数
#ifdef _WINDOWS
CONSOLE_SCREEN_BUFFER_INFO csbi;
GetConsoleScreenBufferInfo(GetStdHandle(STD_OUTPUT_HANDLE), &csbi);
height = csbi.srWindow.Bottom - csbi.srWindow.Top + 1;
#else
struct winsize w;
ioctl(STDOUT_FILENO, TIOCGWINSZ, &w);
height = w.ws_row;
#endif
return height;
}
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;
}
void gctl::parse_filename(std::string filename, std::string &naked_name, std::string &exten_name)
{
naked_name = filename.substr(0, filename.rfind("."));

View File

@ -1,4 +1,4 @@
/********************************************************
/********************************************************
*
*
*
@ -24,104 +24,15 @@
* 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_STREAM_H
#define _GCTL_STREAM_H
// library's head files
#include "stream_t.h"
#ifndef _GCTL_FILE_IO_H
#define _GCTL_FILE_IO_H
// system's head files
#include "cstring"
#include "iostream"
#include "fstream"
#include "iomanip"
#include "ctime"
#include "stdio.h"
#ifdef _WINDOWS
#include "windows.h"
#else
#include "sys/ioctl.h"
#endif
#include "../core/str.h"
namespace gctl
{
/**
* @brief GCTL图标
*/
void display_logo(std::ostream &sout = std::cout);
/**
* @brief
*
* @param new_str
* @param old_str
* @param old_value
* @param new_value
*
* @return
*/
int replace_all(std::string& new_str, const std::string& old_str, const std::string& old_value,
const std::string& new_value);
/**
* @brief
*
* @param in_str
* @param patch_str
* @return
*/
std::string patch_string(std::string in_str, std::string patch_str);
/**
* @brief string对象为stringstream对象
*
* @param[in] in_str string字符串
* @param out_ss stringstreams对象
* @param[in] delimiter string对象时的分隔符
* delimiter转换为空格
*/
void str2ss(std::string in_str, std::stringstream &out_ss, std::string delimiter = "");
/**
* @brief string字符串为double类型的数值
*
* nan或者inf等表示无效值的符号
* >>
* fortran输出文件中以D标注的科学数字转换为浮点类型
*
* @param[in] instr
*
* @return double类型的数值
*/
double str2double(std::string instr);
/**
* @brief double类型数值为string类型字符串 to_string函数
*
* @param[in] in_d
*
* @return
*/
std::string double2str(double in_d);
/**
* @brief
*
* @param[in] length
* @param out
*/
void random_char(unsigned int length, char* out);
/**
* @brief
*
* @param[in] length
* @param out_str
*/
void random_str(unsigned int lenght, std::string &out_str);
/**
/**
* @brief
*
* @param infile ifstream对象
@ -165,30 +76,6 @@ namespace gctl
void open_matched_outfile(std::ofstream &outfile, std::string filename, std::string exten_pattern,
std::ios_base::openmode mode = std::ios_base::out);
/**
* @brief
*
* @return
*/
int terminal_width();
/**
* @brief
*
* @return
*/
int terminal_height();
/**
* @brief
*
* @note
*
* @param in_str
* @param str_vec
*/
void parse_string_with_quotes(std::string in_str, std::vector<std::string> &str_vec);
/**
* @brief
*
@ -197,6 +84,6 @@ namespace gctl
* @param exten_name
*/
void parse_filename(std::string filename, std::string &naked_name, std::string &exten_name);
}
};
#endif // _GCTL_STREAM_H
#endif // _GCTL_FILE_IO_H

View File

@ -29,9 +29,21 @@
#define _GCTL_GMSH_IO_H
// library's head files
#include "../core.h"
#include "../geometry.h"
#include "../utility.h"
#include "../core/matrix.h"
#include "../poly/node.h"
#include "../poly/sphere.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"
@ -101,6 +113,7 @@ 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);
@ -410,6 +423,14 @@ 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)
{
@ -3837,5 +3858,6 @@ namespace gctl
#endif // GCTL_EIGEN
}
};
#endif //_GCTL_GMSH_IO_H

File diff suppressed because it is too large Load Diff

View File

@ -1,526 +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_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 PhysicalTagGeometryTag或者PartitionTag
* @param tag
*/
void edit_group(switch_type_e swt, element_tag_enum tag_type, int tag);
/**
* @brief
*
* @param anchor_type PhysicalTagGeometryTag或者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 PhysicalTagGeometryTag或者PartitionTag
* @param anchor_group
* @param tar_type PhysicalTagGeometryTag或者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 PhysicalTagGeometryTag或者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.2ASCII文件
*
* @param filename
* @param is_packed 0
*/
void read_gmsh_v2_ascii(std::string filename, index_packed_e is_packed = NotPacked);
/**
* @brief Gmsh软件格式的网格剖分文件v2.2ASCII文件
*
* @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

View File

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

View File

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

View File

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

214
lib/io/ply_io.h Normal file
View File

@ -0,0 +1,214 @@
/********************************************************
*
*
*
*
*
*
* 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

325
lib/io/stl_io.h Normal file
View File

@ -0,0 +1,325 @@
/********************************************************
*
*
*
*
*
*
* 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_STL_IO_H
#define _GCTL_STL_IO_H
#include "../poly/triangle.h"
#include "file_io.h"
namespace gctl
{
/**
* @brief Read stl file ascii format
*
* @tparam E vertex attribute type
* @tparam A element attribute type
* @param filename File name
* @param nodes return nodes
* @param facets return facets
* @param solid solid name ('null' will read the first solid)
* @param eps epsilon for node comparison
*/
template <typename E, typename A>
void read_stl_ascii(std::string filename, array<vertex<point3dc, E>> &nodes,
array<type_triangle<A>> &facets, std::string solid = "null", double eps = 1e-8)
{
std::ifstream infile;
open_infile(infile, filename, ".stl");
vertex3dc tmp_node;
std::vector<vertex3dc> tmp_nodes;
_1i_vector one_tri(3);
_2i_vector tmp_tris;
bool inside_solid, new_node;
point3dc nor, vt[3];
char solid_name[256], solid_endname[256];
std::string line;
while (std::getline(infile, line))
{
if (sscanf(line.c_str(), "solid %s", solid_name) == 1 && (!strcmp(solid_name, solid.c_str()) || solid == "null"))
{
inside_solid = true;
while (inside_solid)
{
std::getline(infile, line);
if (sscanf(line.c_str(), "endsolid %s", solid_endname) == 1)
{
if (strcmp(solid_endname, solid_name))
throw std::runtime_error("[gctl::read_stl_ascii] Invalid solid names.");
inside_solid = false;
}
else
{
sscanf(line.c_str(), "facet normal %lf %lf %lf", &nor.x, &nor.y, &nor.z);
std::getline(infile, line); // outer loop
for (int i = 0; i < 3; i++)
{
std::getline(infile, line);
line.erase(0, line.find_first_not_of(" \t"));
line.erase(line.find_last_not_of(" \t") + 1);
sscanf(line.c_str(), "vertex %lf %lf %lf", &vt[i].x, &vt[i].y, &vt[i].z);
}
std::getline(infile, line); // endloop
std::getline(infile, line); // endfacet
for (size_t i = 0; i < 3; i++)
{
tmp_node.id = tmp_nodes.size();
tmp_node.x = vt[i].x;
tmp_node.y = vt[i].y;
tmp_node.z = vt[i].z;
new_node = true;
for (size_t v = 0; v < tmp_nodes.size(); v++)
{
if (isequal(tmp_nodes[v], tmp_node, eps))
{
tmp_node.id = tmp_nodes[v].id;
new_node = false;
break;
}
}
if (new_node) tmp_nodes.push_back(tmp_node);
one_tri[i] = tmp_node.id;
}
tmp_tris.push_back(one_tri);
}
}
}
}
infile.close();
nodes.resize(tmp_nodes.size());
for (size_t i = 0; i < tmp_nodes.size(); i++)
{
nodes[i].id = tmp_nodes[i].id;
nodes[i].x = tmp_nodes[i].x;
nodes[i].y = tmp_nodes[i].y;
nodes[i].z = tmp_nodes[i].z;
}
facets.resize(tmp_tris.size());
for (size_t i = 0; i < tmp_tris.size(); i++)
{
facets[i].set(nodes[tmp_tris[i][0]], nodes[tmp_tris[i][1]], nodes[tmp_tris[i][2]], i);
}
destroy_vector(tmp_nodes);
destroy_vector(tmp_tris);
destroy_vector(one_tri);
return;
}
/**
* @brief Save stl file ascii format
*
* @tparam E Vertex attribute type
* @tparam A Element attribute type
* @param filename File name
* @param facets Input facets
* @param solid Solid name
*/
template <typename A>
void save_stl_ascii(std::string filename, const array<type_triangle<A>> &facets,
std::string solid = "STL generated by GCTL")
{
std::ofstream outfile;
open_outfile(outfile, filename, ".stl");
point3dc nor;
outfile << "solid " << solid << "\n";
for (size_t i = 0; i < facets.size(); i++)
{
nor = cross(*facets[i].vert[1] - *facets[i].vert[0], *facets[i].vert[2] - *facets[i].vert[0]).normal();
outfile << "facet normal " << nor.x << " " << nor.y << " " << nor.z << "\n";
outfile << " outer loop\n";
outfile << " vertex " << facets[i].vert[0]->x << " " << facets[i].vert[0]->y << " " << facets[i].vert[0]->z << "\n";
outfile << " vertex " << facets[i].vert[1]->x << " " << facets[i].vert[1]->y << " " << facets[i].vert[1]->z << "\n";
outfile << " vertex " << facets[i].vert[2]->x << " " << facets[i].vert[2]->y << " " << facets[i].vert[2]->z << "\n";
outfile << " endloop\n";
outfile << "endfacet\n";
}
outfile << "endsolid " << solid << "\n";
outfile.close();
return;
}
/**
* @brief Read stl file binary format
*
* @tparam E vertex attribute type
* @tparam A element attribute type
* @param filename File name
* @param nodes return nodes
* @param facets return facets
* @param eps epsilon for node comparison
*
* @return solid name
*/
template <typename E, typename A>
std::string read_stl_binary(std::string filename, array<vertex<point3dc, E>> &nodes,
array<type_triangle<A>> &facets, double eps = 1e-8)
{
std::ifstream infile;
open_infile(infile, filename, ".stl", std::ios::binary);
char blank[2], header[80];
infile.read(header, 80);
int n_facets;
infile.read((char*)&n_facets, 4);
vertex3dc tmp_node;
std::vector<vertex3dc> tmp_nodes;
_1i_vector one_tri(3);
_2i_vector tmp_tris;
tmp_tris.reserve(n_facets);
bool new_node;
point3fc nor, vt[3];
for (size_t i = 0; i < n_facets; i++)
{
infile.read((char*)&nor.x, 4);
infile.read((char*)&nor.y, 4);
infile.read((char*)&nor.z, 4);
for (size_t j = 0; j < 3; j++)
{
infile.read((char*)&vt[j].x, 4);
infile.read((char*)&vt[j].y, 4);
infile.read((char*)&vt[j].z, 4);
}
infile.read((char*)&blank, 2);
for (size_t j = 0; j < 3; j++)
{
tmp_node.id = tmp_nodes.size();
tmp_node.x = vt[j].x;
tmp_node.y = vt[j].y;
tmp_node.z = vt[j].z;
new_node = true;
for (size_t v = 0; v < tmp_nodes.size(); v++)
{
if (isequal(tmp_nodes[v], tmp_node, eps))
{
tmp_node.id = tmp_nodes[v].id;
new_node = false;
break;
}
}
if (new_node) tmp_nodes.push_back(tmp_node);
one_tri[j] = tmp_node.id;
}
tmp_tris.push_back(one_tri);
}
infile.close();
nodes.resize(tmp_nodes.size());
for (size_t i = 0; i < tmp_nodes.size(); i++)
{
nodes[i].id = tmp_nodes[i].id;
nodes[i].x = tmp_nodes[i].x;
nodes[i].y = tmp_nodes[i].y;
nodes[i].z = tmp_nodes[i].z;
}
facets.resize(tmp_tris.size());
for (size_t i = 0; i < tmp_tris.size(); i++)
{
facets[i].set(nodes[tmp_tris[i][0]], nodes[tmp_tris[i][1]], nodes[tmp_tris[i][2]], i);
}
destroy_vector(tmp_nodes);
destroy_vector(tmp_tris);
destroy_vector(one_tri);
std::string solid_name(header);
return solid_name;
}
/**
* @brief Save stl file binary format
*
* @tparam E vertex attribute type
* @tparam A element attribute type
* @param filename File name
* @param facets return facets
* @param solid solid name
*/
template <typename A>
void save_stl_binary(std::string filename, const array<type_triangle<A>> &facets,
std::string solid = "STL generated by GCTL")
{
std::ofstream outfile;
open_outfile(outfile, filename, ".stl", std::ios::binary);
char blank[2] = {'\0', '\0'};
char header[80];
for (size_t i = 0; i < 80; i++) header[i] = '\0';
for (size_t i = 0; i < solid.size(); i++) header[i] = solid[i];
outfile.write(header, 80);
int n_facets = facets.size();
outfile.write((char*)&n_facets, 4);
point3dc nor;
point3fc nor_fl;
point3fc vt[3];
for (size_t i = 0; i < facets.size(); i++)
{
nor = cross(*facets[i].vert[1] - *facets[i].vert[0], *facets[i].vert[2] - *facets[i].vert[0]).normal();
// convert from double to float
nor_fl.x = nor.x;
nor_fl.y = nor.y;
nor_fl.z = nor.z;
for (size_t j = 0; j < 3; j++)
{
vt[j].x = facets[i].vert[j]->x;
vt[j].y = facets[i].vert[j]->y;
vt[j].z = facets[i].vert[j]->z;
}
outfile.write((char*)&nor_fl.x, 4);
outfile.write((char*)&nor_fl.y, 4);
outfile.write((char*)&nor_fl.z, 4);
for (size_t j = 0; j < 3; j++)
{
outfile.write((char*)&vt[j].x, 4);
outfile.write((char*)&vt[j].y, 4);
outfile.write((char*)&vt[j].z, 4);
}
outfile.write(blank, 2);
}
outfile.close();
}
};
#endif // _GCTL_STL_IO_H

View File

@ -29,8 +29,8 @@
#define _GCTL_SURFER_IO_H
// library's head files
#include "../core.h"
#include "../utility.h"
#include "../core/matrix.h"
#include "file_io.h"
namespace gctl
{
@ -187,6 +187,6 @@ namespace gctl
void save_surfer7_grid(std::string filename, const array<double> &out_data,
int xnum, int ynum, double xmin, double ymin, double dx, double dy,
double zmin = NAN, double zmax = NAN, double blank_val = 1.70141e38);
}
};
#endif // _GCTL_SURFER_IO_H

170
lib/io/term_io.cpp Normal file
View File

@ -0,0 +1,170 @@
/********************************************************
*
*
*
*
*
*
* 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.
******************************************************/
#if defined _WINDOWS || __WIN32__
#include <io.h>
// Test for file existence
#define F_OK 0
#endif
#include "term_io.h"
//在终端显示一个简易的GCTL图标
void gctl::display_logo(std::ostream &sout)
{
sout << " ___ ___ _____ __\n";
sout << " / _ \\ / __\\/__ \\ / /\n";
sout << " / /_\\// / / /\\// /\n";
sout << "/ /_\\\\/ /___ / / / /___\n";
sout << "\\____/\\____/ \\/ \\____/\n";
sout << "Geophysical Computational Tools & Library\n";
sout << "Author: Dr. Yi Zhang (yizhang-geo@zju.edu.cn)\n\n";
return;
}
int gctl::terminal_width()
{
int width;
//获取终端窗口的行列数
#ifdef _WINDOWS
CONSOLE_SCREEN_BUFFER_INFO csbi;
GetConsoleScreenBufferInfo(GetStdHandle(STD_OUTPUT_HANDLE), &csbi);
width = csbi.srWindow.Right - csbi.srWindow.Left + 1;
#else
struct winsize w;
ioctl(STDOUT_FILENO, TIOCGWINSZ, &w);
width = w.ws_col;
#endif
return width;
}
int gctl::terminal_height()
{
int height;
//获取终端窗口的行列数
#ifdef _WINDOWS
CONSOLE_SCREEN_BUFFER_INFO csbi;
GetConsoleScreenBufferInfo(GetStdHandle(STD_OUTPUT_HANDLE), &csbi);
height = csbi.srWindow.Bottom - csbi.srWindow.Top + 1;
#else
struct winsize w;
ioctl(STDOUT_FILENO, TIOCGWINSZ, &w);
height = w.ws_row;
#endif
return height;
}
gctl::cli_viewer::cli_viewer()
{
// 初始化ncurses
initscr();
cbreak(); // 禁用行缓冲
noecho(); // 不显示输入字符
keypad(stdscr, TRUE); // 启用键盘扩展键
timeout(0); // 设置非阻塞输入
curs_set(0); // 隐藏光标
}
// 析构函数
gctl::cli_viewer::~cli_viewer()
{
// 清理ncurses
endwin();
}
// 设置显示内容
void gctl::cli_viewer::setData(const std::vector<std::string>& data)
{
lines = data;
return;
}
void gctl::cli_viewer::addData(const std::string& l)
{
lines.push_back(l);
return;
}
// 显示文件内容并进入查看循环
void gctl::cli_viewer::display()
{
int maxLines = LINES - 1; // 终端的行数(减去提示行)
int startLine = 0; // 当前显示的起始行
int startX = 0; // 当前水平滚动的起始位置
int maxLineLength = 0; // 最长行的长度
// 计算最长行的长度
for (const auto& line : lines)
{
maxLineLength = std::max(maxLineLength, static_cast<int>(line.length()));
}
// 主循环
while (true)
{
werase(stdscr); // 清屏但不闪烁
// 显示文件内容
for (int i = startLine; i < startLine + maxLines && i < lines.size(); ++i)
{
const std::string& line = lines[i];
if (startX < line.length())
{
mvprintw(i - startLine, 0, "%s", line.substr(startX).c_str());
}
}
// 在最后一行显示提示信息
mvprintw(LINES - 1, 0, "Use arrow keys or mouse wheel to scroll. Press Q or ESC to exit.");
int ch = getch(); // 获取按键
if (ch == KEY_UP)
{
startLine = std::max(0, startLine - 1); // 向上滚动
}
else if (ch == KEY_DOWN)
{
startLine = std::min(startLine + 1, static_cast<int>(lines.size()) - maxLines); // 向下滚动
}
else if (ch == KEY_LEFT)
{
startX = std::max(0, startX - 1); // 向左滚动
}
else if (ch == KEY_RIGHT)
{
startX = std::min(startX + 1, maxLineLength - COLS); // 向右滚动
}
else if (ch == 'q' || ch == 'Q' || ch == 27)
{ // 按Q或ESC退出
break;
}
wrefresh(stdscr); // 刷新屏幕
}
return;
}

103
lib/io/term_io.h Normal file
View File

@ -0,0 +1,103 @@
/********************************************************
*
*
*
*
*
*
* 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_TERM_IO_H
#define _GCTL_TERM_IO_H
#include <iostream>
#include <fstream>
#include <vector>
#include <unistd.h>
#include <stdio.h>
#include <ncurses.h>
#ifdef _WINDOWS
#include <windows.h>
#else
#include <sys/ioctl.h>
#endif
/**
* terminal control symbols
*/
#define GCTL_BOLD "\033[1m" ///< 设置终端字体与颜色为粗体
#define GCTL_BOLDRED "\033[1m\033[31m" ///< 设置终端字体与颜色为粗体红色
#define GCTL_BOLDGREEN "\033[1m\033[32m" ///< 设置终端字体与颜色为粗体绿色
#define GCTL_BOLDBLUE "\033[1m\033[34m" ///< 设置终端字体与颜色为粗体蓝色
#define GCTL_BOLDYELLOW "\033[1m\033[33m" ///< 设置后续字符字体为黄色加粗
#define GCTL_UNDERLINE "\033[1m\033[4m" ///< 设置终端字体为下划线
#define GCTL_RESET "\033[0m" ///< 重置字体与颜色设定
#define GCTL_CLEARLINE "\033[K" ///< 清空终端当前行
#define GCTL_CLEARALL "\033[2J" ///< 清空终端窗口
#define GCTL_MOVEUP(os, x) do {os << "\033[" << x << "A";} while (0); ///< 终端光标上移x行
#define GCTL_MOVEDOWN(os, x) do {os << "\033[" << x << "B";} while (0); ///< 终端光标下移x行
#define GCTL_MOVERIGHT(os, x) do {os << "\033[" << x << "C";} while (0); ///< 终端光标右移x列
#define GCTL_MOVELEFT(os, x) do {os << "\033[" << x << "D";} while (0); ///< 终端光标左移x列
#define GCTL_MOVETO(os, r, c) do {os << "\033[" << c << ";" << r << "H";} while (0); ///< 终端光标移动至r行c列位置
namespace gctl
{
/**
* @brief GCTL图标
*/
void display_logo(std::ostream &sout = std::cout);
/**
* @brief
*
* @return
*/
int terminal_width();
/**
* @brief
*
* @return
*/
int terminal_height();
// 命令行查看器,用于在终端内查看文本文件内容,支持上下左右鼠标滚动或按键导航和退出查看
class cli_viewer
{
public:
// 构造函数
cli_viewer();
// 析构函数
~cli_viewer();
// 设置显示内容
void setData(const std::vector<std::string>& data);
// 添加显示内容
void addData(const std::string& l);
// 显示文件内容并进入查看循环
void display();
private:
std::vector<std::string> lines; // 存储文件内容
};
};
#endif // _GCTL_TERM_IO_H

View File

@ -28,10 +28,11 @@
#ifndef _GCTL_TETGEN_IO_H
#define _GCTL_TETGEN_IO_H
// library's head file
#include "../core.h"
#include "../utility.h"
#include "../geometry.h"
#include "../core/matrix.h"
#include "../poly/triangle.h"
#include "../poly/tetrahedron.h"
#include "../poly/edge.h"
#include "file_io.h"
namespace gctl
{
@ -1192,5 +1193,6 @@ namespace gctl
tetout.close();
return;
}
}
};
#endif //_GCTL_TETGEN_IO_H

View File

@ -28,11 +28,11 @@
#ifndef _GCTL_TEXT_IO_H
#define _GCTL_TEXT_IO_H
// library's head file
#include "../core.h"
#include "../utility.h"
#include "../geometry.h"
// system head file
#include "../core/matrix.h"
#include "../poly/point3c.h"
#include "../poly/point3s.h"
#include "file_io.h"
#include "initializer_list"
namespace gctl
@ -1322,6 +1322,7 @@ namespace gctl
*/
void save_data_column(std::ofstream &outfile, std::initializer_list<std::vector<double>*> dat_val,
std::initializer_list<std::string> dat_name, char delimiter = ' ', int precision = GCTL_PRECISION);
}
};
#endif //_GCTL_TEXT_IO_H

View File

@ -29,9 +29,11 @@
#define _GCTL_TRIANGLE_IO_H
// library's head file
#include "../core.h"
#include "../utility.h"
#include "../geometry.h"
#include "../core/matrix.h"
#include "../poly/triangle2d.h"
#include "../poly/triangle2d2o.h"
#include "../poly/edge2d.h"
#include "file_io.h"
namespace gctl
{
@ -1270,5 +1272,7 @@ namespace gctl
triout.close();
return;
}
}
};
#endif //_GCTL_TRIANGLE_H

View File

@ -25,40 +25,43 @@
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#include "refellipsoid.h"
#ifndef _GCTL_MATH_H
#define _GCTL_MATH_H
gctl::refellipsoid::refellipsoid(/* args */){}
#include "math/gmath.h"
gctl::refellipsoid::refellipsoid(refellipsoid_type_e refellipsoid)
{
set(refellipsoid);
}
#include "math/legendre.h"
#include "math/linear_algebra.h"
gctl::refellipsoid::refellipsoid(double R, double r, bool flat)
{
set(R, r, flat);
}
#include "math/heap_sort.h"
#include "math/boxsort2d.h"
#include "math/boxsort_sph.h"
#include "math/variogram.h"
gctl::refellipsoid::~refellipsoid(){}
#include "math/interpolate.h"
#include "math/extrapolate.h"
void gctl::refellipsoid::set(refellipsoid_type_e refellipsoid)
{
if (refellipsoid == Earth) {r_ = R_ = GCTL_Earth_Radius;}
else if (refellipsoid == Moon) {r_ = R_ = GCTL_Moon_Radius;}
else if (refellipsoid == Mars) {r_ = R_ = GCTL_Mars_Radius;}
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/gaussfunc.h"
#include "math/sinkhorn.h"
#include "math/kde.h"
#include "math/autodiff.h"
#include "math/multinary.h"
#include "math/emd.h"
void gctl::refellipsoid::set(double R, double r_or_flat, bool is_flat)
{
R_ = R;
if (is_flat) r_ = R_*(1 - 1/r_or_flat);
else r_ = r_or_flat;
}
#include "math/geometry2d.h"
#include "math/geometry3d.h"
#include "math/refellipsoid.h"
double gctl::refellipsoid::radius_at(double lati)
{
return ellipse_radius_2d(R_, r_, lati*M_PI/180.0, 0.0);
}
#include "math/mathfunc_ext.h"
#include "math/space_filter.h"
#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

View File

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

View File

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

View File

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

View File

@ -29,8 +29,6 @@
#define _GCTL_EMD_H
// library head file
#include "../gctl_config.h"
#include "../core/array.h"
#include "../core/matrix.h"
#ifdef GCTL_EEMD
@ -46,7 +44,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

View File

@ -28,9 +28,7 @@
#ifndef _GCTL_EXTRAPOLATE_H
#define _GCTL_EXTRAPOLATE_H
#include "../core/array.h"
#include "../core/matrix.h"
#include "../maths/mathfunc_t.h"
#include "gmath.h"
namespace gctl
{
@ -105,6 +103,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

View File

@ -54,6 +54,7 @@ void gctl::dft_r2c_1d(const _1d_array &in_real, _1cd_array &out_spectrum, double
out_spectrum.resize(n);
for (int i = 0; i < n; i++)
{
// 这里我们把2.0改成了1.0
out_spectrum[i].real(out[i][0]*2.0/m);
out_spectrum[i].imag(out[i][1]*2.0/m);
}
@ -88,6 +89,7 @@ void gctl::dft_c2r_1d(const _1cd_array &in_spectrum, _1d_array &out_real)
for (int i = 0; i < m; i++)
{
// 这里我们删除了 *0.5
in[i][0] = in_spectrum[i].real()*0.5*n;
in[i][1] = in_spectrum[i].imag()*0.5*n;
}
@ -223,4 +225,110 @@ void gctl::idct1d(const _1d_array &in_arr, _1d_array &out_arr)
return;
}
#else
void gctl::dft_r2c_1d(const _1d_array &in_real, _1cd_array &out_spectrum, double sampling, _1d_array *freq_ptr)
{
if (in_real.empty())
{
throw std::invalid_argument("[GCTL] Invalid input array size for gctl::dft_r2c_1d(...)");
}
int n = in_real.size();
int s = n/2 + 1; // out size
int m = log2(n);
double *ar, *ai;
ar = new double [n];
ai = new double [n];
for (size_t i = 0; i < n; i++)
{
ar[i] = in_real[i];
ai[i] = 0.0;
}
int i,j,k,k1,l1,l2,ip;
double *ir,*ii,*i1,*i2;
double f1,t1,t2,v,u,w1,w2,d,vu,pi1;
f1=(double)1;
j=0;
for(i=0,ir=ar,ii=ai;i<n;i++,ir++,ii++)
{ if(i<j)
{ i1=ar+j; i2=ai+j;
t1=*i1; t2=*i2; *i1=*ir; *i2=*ii; *ir=t1; *ii=t2;
}
k1=n/2;
while((j>=k1)&&(k1!=0))
{ j=j-k1; k1=k1/2;
}
j=j+k1;
}
for(j=1;j<m+1;j++)
{ l1=2;
for (i=1;i<j;i++) l1=2*l1;
l2=l1/2; pi1=(double)GCTL_Pi/(double)l2; v=(double)1; u=(double)0;
w1=(double)cos(pi1); w2=(double)sin(pi1);
for(k=1;k<l2+1;k++)
{ d=f1*u;
for(i=k-1;i<n-l1+k;i+=l1)
{ ip=i+l2; i1=ar+ip; i2=ai+ip; ir=ar+i; ii=ai+i;
t1=*i1*v+*i2*d; t2=*i2*v-*i1*d;
*i1=*ir-t1; *i2=*ii-t2; *ir=*ir+t1; *ii=*ii+t2;
}
vu=v*w1-u*w2; u=v*w2+u*w1; v=vu;
} }
for(i=0,ir=ar,ii=ai;i<n;i++,ir++,ii++)
{
*ir/=n; *ii/=n;
}
out_spectrum.resize(s);
for (int i = 0; i < s; i++)
{
out_spectrum[i].real(ar[i]);
out_spectrum[i].imag(ai[i]);
}
if (freq_ptr != nullptr)
{
freq_ptr->resize(s);
for (size_t i = 0; i < s; i++)
{
freq_ptr->at(i) = 0.5*sampling*i/s;
}
}
delete[] ar;
delete[] ai;
return;
}
void gctl::dft_c2r_1d(const _1cd_array &in_spectrum, _1d_array &out_real)
{
return;
}
void gctl::dft2d(const _2d_matrix &in_arr, _1cd_array &out_arr)
{
return;
}
void gctl::idft2d(const _1cd_array &in_arr, int m, int n, _2d_matrix &out_arr)
{
return;
}
void gctl::dct1d(const _1d_array &in_arr, _1d_array &out_arr)
{
return;
}
void gctl::idct1d(const _1d_array &in_arr, _1d_array &out_arr)
{
return;
}
#endif // GCTL_FFTW3

View File

@ -28,16 +28,11 @@
#ifndef _GCTL_FFT_H
#define _GCTL_FFT_H
// library head file
#include "../gctl_config.h"
#include "../core/array.h"
#include "../core/matrix.h"
#include "../core/exceptions.h"
#ifdef GCTL_FFTW3
// fftw head file
#include "fftw3.h"
#endif // GCTL_FFTW3
namespace gctl
{
@ -92,8 +87,6 @@ namespace gctl
* @param out_arr Output real array
*/
void idct1d(const _1d_array &in_arr, _1d_array &out_arr);
}
#endif // GCTL_FFTW3
};
#endif // _GCTL_FFT_H

View File

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

View File

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

View File

@ -77,6 +77,6 @@ namespace gctl
*/
double gaussian_dist2d(double x, double y, const gaussian_para2d &p,
gradient_type_e mode = Value);
}
};
#endif // _GAUSSFUNC_H

View File

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

View File

@ -28,14 +28,8 @@
#ifndef _GCTL_GEOMETRY2D_H
#define _GCTL_GEOMETRY2D_H
// library's head files
#include "../core/array.h"
#include "../maths.h"
#include "triangle2d.h"
#include "triangle2d2o.h"
#include "edge2d.h"
#include "rectangle2d.h"
#include "../poly/triangle2d.h"
#include "../poly/edge2d.h"
namespace gctl
{
@ -239,5 +233,6 @@ 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

View File

@ -25,15 +25,16 @@
* 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()));
@ -110,7 +111,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 = p1.s2c(), c2 = p2.s2c();
gctl::point3dc c1 = s2c(p1), c2 = s2c(p2);
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;
}
@ -135,11 +136,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 = v1.s2c(), p2 = v2.s2c();
point3dc p1 = s2c(v1), p2 = s2c(v2);
double Phi = angle(p1, p2); // Phi为弧度表示
point3dc pc = cos(arc)*p1 + sin(arc)/sin(Phi)*(p2 - cos(Phi)*p1);
return pc.c2s();
return c2s(pc);
}
gctl::point3dc gctl::geometry3d::dot_on_plane(const point3dc &face_dot, const point3dc &face_nor,

View File

@ -28,21 +28,9 @@
#ifndef _GCTL_GEOMETRY3D_H
#define _GCTL_GEOMETRY3D_H
// library's head file
#include "../core/array.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"
#include "../poly/triangle.h"
#include "../poly/tetrahedron.h"
#include "gmath.h"
namespace gctl
{
@ -361,5 +349,6 @@ 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

View File

@ -66,7 +66,7 @@ void glni::coefficients::initiate(size_t size)
#ifdef GCTL_CHECK_SIZE
if (size == 0 || size > 20)
{
throw invalid_argument("Invalid initiating order. From glni::coefficients::initiate(...).");
throw std::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 out_of_range("Invalid index. From gctl::glni::coefficients::node(...)");
throw std::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 out_of_range("Invalid index. From gctl::glni::coefficients::weight(...)");
throw std::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 invalid_argument("Invalid initiating order. From glni::line::initiate(...)");
throw std::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 length_error("Invalid integration range. From glni::line::integral(...)");
throw std::length_error("Invalid integration range. From glni::line::integral(...)");
}
#endif
if (!initialized_)
{
throw runtime_error("The object is not initialized. From gctl::glni::line::integral(...).");
throw std::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 invalid_argument("Invalid initiating order. From glni::quadrilateral::initiate(...)");
throw std::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 length_error("Invalid integration range. From glni::quadrilateral::integral(...)");
throw std::length_error("Invalid integration range. From glni::quadrilateral::integral(...)");
}
#endif
if (!initialized_)
{
throw runtime_error("Un-initialized. From glni::quadrilateral::integral(...)");
throw std::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 invalid_argument("Invalid initiating order. From glni::triangle::initiate(...)");
throw std::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 length_error("Invalid integration range. From glni::triangle::integral(...)");
throw std::length_error("Invalid integration range. From glni::triangle::integral(...)");
}
#endif
if (!initialized_)
{
throw runtime_error("Un-initialized. From glni::triangle::integral(...)");
throw std::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 invalid_argument("Invalid initiating order. From glni::triangle_c::initiate(...)");
throw std::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) <= GCTL_ZERO) // 这里没有设置为0是为了检测三个点是否共线
if ((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1) <= 1e-30) // 这里没有设置为0是为了检测三个点是否共线
{
throw length_error("Invalid integration range. From glni::triangle::integral(...)");
throw std::length_error("Invalid integration range. From glni::triangle::integral(...)");
}
#endif
if (!initialized_)
{
throw runtime_error("Un-initialized. From glni::triangle::integral(...)");
throw std::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 invalid_argument("Invalid initiating order. From glni::cube::initiate(...)");
throw std::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 length_error("Invalid integration range. From glni::cube::integral(...)");
throw std::length_error("Invalid integration range. From glni::cube::integral(...)");
}
#endif
if (!initialized_)
{
throw runtime_error("Un-initialized. From glni::cube::integral(...)");
throw std::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 invalid_argument("Invalid initiating order. From glni::tetrahedron::initiate(...)");
throw std::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 runtime_error("Un-initialized. From glni::tetrahedron::integral(...)");
throw std::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 < GCTL_ZERO)
if (scl < 1e-30)
{
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::invalid_argument("glni::triprism::initiate(...) got an invalid initiating order.");
throw std::std::invalid_argument("glni::triprism::initiate(...) got an invalid initiating order.");
}
#endif // GCTL_CHECK_BOUNDER

View File

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

View File

@ -1,4 +1,4 @@
/********************************************************
/********************************************************
*
*
*
@ -25,7 +25,483 @@
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#include "algorithm_func.h"
#include "gmath.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)
{

View File

@ -25,21 +25,145 @@
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#ifndef _GCTL_MATHFUNC_H
#define _GCTL_MATHFUNC_H
#ifndef _GCTL_GMATH_H
#define _GCTL_GMATH_H
#include <cmath>
#include <random>
#include <vector>
// library's head files
#include "../core/macro.h"
#include "../core/spmat.h"
#include "mathfunc_t.h"
// system's head files
#include "cmath"
#include "random"
#include "vector"
namespace gctl
{
template <typename T>
inline int sign(T d)
{
return (T(0) < d) - (d < T(0));
}
template <typename T>
inline T arc(T deg)
{
return deg*GCTL_Pi/180.0L;
}
template <typename T>
inline T deg(T arc)
{
return arc*180.0L/GCTL_Pi;
}
template <typename T>
inline T sind(T deg)
{
return sin(deg*GCTL_Pi/180.0L);
}
template <typename T>
inline T cosd(T deg)
{
return cos(deg*GCTL_Pi/180.0L);
}
template <typename T>
inline T tand(T deg)
{
return tan(deg*GCTL_Pi/180.0L);
}
template <typename T>
inline T power2(T in)
{
return (in)*(in);
}
template <typename T>
inline T power3(T in)
{
return (in)*(in)*(in);
}
template <typename T>
inline T power4(T in)
{
return (in)*(in)*(in)*(in);
}
template <typename T>
inline T power5(T in)
{
return (in)*(in)*(in)*(in)*(in);
}
template <typename T>
inline T jacoby2(T x00, T x01, T x10, T x11)
{
return x00*x11-x01*x10;
}
template <typename T>
inline T jacoby3(T x00, T x01, T x02, T x10, T x11,
T x12, T x20, T x21, T x22)
{
return x00*x11*x22+x01*x12*x20+x10*x21*x02
-x02*x11*x20-x01*x10*x22-x00*x12*x21;
}
/**
* @brief Calculate the inverse matrix of a 3x3 matrix
*
* @tparam T Type name
* @param A Pointer of the input matrix, must be stored in an array of the nine coefficients in a row-major fashion
* @param invA Pointer of the output matrix, must be stored in an array of the nine coefficients in a row-major fashion
* @return true Success
* @return false Fail
*/
template <typename T>
bool inverse3x3(T *A, T *invA)
{
T det = jacoby3(A[0], A[1], A[2], A[3], A[4], A[5], A[6], A[7], A[8]);
if (det <= GCTL_ZERO && det >= -1*GCTL_ZERO) return false;
invA[0] = jacoby2(A[4], A[7], A[5], A[8])/det;
invA[1] = -1.0*jacoby2(A[1], A[7], A[2], A[8])/det;
invA[2] = jacoby2(A[1], A[4], A[2], A[5])/det;
invA[3] = -1.0*jacoby2(A[3], A[6], A[5], A[8])/det;
invA[4] = jacoby2(A[0], A[6], A[2], A[8])/det;
invA[5] = -1.0*jacoby2(A[0], A[3], A[2], A[5])/det;
invA[6] = jacoby2(A[3], A[6], A[4], A[7])/det;
invA[7] = -1.0*jacoby2(A[0], A[6], A[1], A[7])/det;
invA[8] = jacoby2(A[0], A[3], A[1], A[4])/det;
return true;
}
template <typename T>
T arctg(T v)
{
T ang;
if(v>=0) ang=atan(v);
else ang=atan(v)+GCTL_Pi;
return ang;
}
template <typename T>
T arctg2(T v, T f)
{
T ang;
if(f>=0)
{
if(atan(v)>0) ang=atan(v);
else ang=atan(v) + GCTL_Pi;
}
else if(f<0)
{
if(atan(v)<0) ang=atan(v);
else ang=atan(v) - GCTL_Pi;
}
return ang;
}
/**
* @brief Return a random number in the range [low, hig]
*
@ -212,6 +336,107 @@ namespace gctl
* @param[in] a_s
*/
void schmidt_orthogonal(const array<double> &a, array<double> &e, int a_s);
}
#endif // _GCTL_MATHFUNC_H
/**
* @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 14使
*/
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 14使
*/
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 14使
*/
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_GMATH_H

View File

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

View File

@ -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 = out_posi[i].s2c();
tmp_pc2 = in_posi[tmp_id].s2c();
tmp_pc1 = s2c(out_posi[i]);
tmp_pc2 = s2c(in_posi[tmp_id]);
tmp_dist = distance(tmp_pc1, tmp_pc2);
if (geometry3d::angle(tmp_pc1, tmp_pc2) <= search_deg &&

View File

@ -28,13 +28,7 @@
#ifndef _GCTL_INTERPOLATE_H
#define _GCTL_INTERPOLATE_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"
#include "geometry3d.h"
namespace gctl
{
@ -316,6 +310,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

View File

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

View File

@ -28,7 +28,7 @@
#ifndef _GCTL_LEGENDRE_H
#define _GCTL_LEGENDRE_H
#include "../core.h"
#include "../core/array.h"
namespace gctl
{
@ -98,6 +98,6 @@ namespace gctl
void nalf_sfcm(array<array<double>> &nalf, const array<array<double>> &a_nm,
const array<array<double>> &b_nm, int max_order, double theta,
legendre_norm_e norm, bool derivative = false);
}
};
#endif //_GCTL_LEGENDRE_H

View File

@ -28,8 +28,6 @@
#ifndef _GCTL_LINEAR_ALGEBRA_H
#define _GCTL_LINEAR_ALGEBRA_H
#include "../core.h"
#ifdef GCTL_OPENMP
#include "omp.h"
#endif // GCTL_OPENMP
@ -38,6 +36,8 @@
#include "cblas.h"
#endif // GCTL_OPENBLAS
#include "../core/matrix.h"
namespace gctl
{
/**
@ -350,6 +350,6 @@ namespace gctl
* @return New standard deviation value.
*/
double dynamic_stddev(double old_stddev, size_t old_count, double old_mean, double new_input, double &new_mean);
}
};
#endif // _GCTL_LINEAR_ALGEBRA_H

View File

@ -30,6 +30,10 @@
namespace gctl
{
/**
*
*/
void DTfunc(double **a,int m,int n,double A,double B,int sw);
void DTari(double **a,int m,int n,double A,int sw);
void DTfilter(double **a,int m,int n,double A,double B,double Q,int sw);
@ -47,6 +51,6 @@ namespace gctl
void tran(double **a,int n);
void metrimult(double **a,double **b,int m,int n,int l,double **c);
void svd_standalone(double **a,int m,int n,double **u,double **v,double *q,int Index);
}
};
#endif // _GCTL_MATHFUNC_EXT_H

View File

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

143
lib/math/refellipsoid.cpp Normal file
View File

@ -0,0 +1,143 @@
/********************************************************
*
*
*
*
*
*
* 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 "refellipsoid.h"
gctl::refellipsoid::refellipsoid(/* args */){}
gctl::refellipsoid::refellipsoid(refellipsoid_type_e refellipsoid)
{
set(refellipsoid);
}
gctl::refellipsoid::refellipsoid(double R, double r, bool flat)
{
set(R, r, flat);
}
gctl::refellipsoid::~refellipsoid(){}
void gctl::refellipsoid::set(refellipsoid_type_e refellipsoid)
{
if (refellipsoid == Earth) {r_ = R_ = GCTL_Earth_Radius;}
else if (refellipsoid == MagEarth) {r_ = R_ = GCTL_Earth_RefRadius;}
else if (refellipsoid == Moon) {r_ = R_ = GCTL_Moon_Radius;}
else if (refellipsoid == Mars) {r_ = R_ = GCTL_Mars_Radius;}
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(...)");
f_ = (R_ - r_)/R_;
e_ = sqrt(1.0 - (r_*r_)/(R_*R_));
return;
}
void gctl::refellipsoid::set(double R, double r_or_flat, bool is_flat)
{
R_ = R;
if (is_flat) r_ = R_*(1.0 - 1.0/r_or_flat);
else r_ = r_or_flat;
f_ = (R_ - r_)/R_;
e_ = sqrt(1.0 - (r_*r_)/(R_*R_));
return;
}
double gctl::refellipsoid::geodetic_radius(double geodetic_lati)
{
return R_/sqrt(1.0 - e_*e_*sind(geodetic_lati)*sind(geodetic_lati));
}
void gctl::refellipsoid::geodetic2spherical(double geodetic_lati,
double geodetic_hei, double& sph_lati, double& sph_rad)
{
double CosLat, SinLat, rc, xp, zp;
/*
** Convert geodetic coordinates, (defined by the WGS-84
** reference ellipsoid), to Earth Centered Earth Fixed Cartesian
** coordinates, and then to spherical coordinates.
*/
CosLat = cosd(geodetic_lati);
SinLat = sind(geodetic_lati);
// compute the local rRdius of curvature on the WGS-84 reference ellipsoid
rc = R_/sqrt(1.0 - e_*e_*SinLat*SinLat);
// compute ECEF Cartesian coordinates of specified point (for longitude=0)
xp = (rc + geodetic_hei)*CosLat;
zp = (rc*(1.0 - e_*e_) + geodetic_hei)*SinLat;
// compute spherical rRdius and angle lambda and phi of specified point
sph_rad = sqrt(xp*xp + zp*zp);
sph_lati = deg(asin(zp/sph_rad));
return;
}
void gctl::refellipsoid::spherical2geodetic(const point3ds& ps, double& geodetic_lon, double& geodetic_lati,
double& geodetic_hei, double eps, int cnt)
{
point3dc pc = s2c(ps);
xyz2geodetic(pc, geodetic_lon, geodetic_lati, geodetic_hei, eps, cnt);
return;
}
void gctl::refellipsoid::xyz2geodetic(const point3dc& pc, double& geodetic_lon,
double& geodetic_lati, double& geodetic_hei, double eps, int cnt)
{
double curB, N;
double H = sqrt(pc.x*pc.x + pc.y*pc.y);
double calB = atan2(pc.z, H);
int c = 0;
do
{
curB = calB;
N = R_/sqrt(1 - e_*e_*sin(curB)*sin(curB));
calB = atan2(pc.z + N*e_*e_*sin(curB), H);
c++;
}
while (abs(curB - calB)*180.0/GCTL_Pi > eps && c < cnt);
geodetic_lon = atan2(pc.y, pc.x)*180.0/GCTL_Pi;
geodetic_lati = curB*180.0/GCTL_Pi;
geodetic_hei = pc.z/sin(curB) - N*(1 - e_*e_);
return;
}
void gctl::refellipsoid::geodetic2xyz(double geodetic_lon, double geodetic_lati,
double geodetic_hei, point3dc& pc)
{
double L = arc(geodetic_lon);
double B = arc(geodetic_lati);
double H = geodetic_hei;
double N = R_/sqrt(1 - e_*e_*sin(B)*sin(B));
pc.x = (N + H)*cos(B)*cos(L);
pc.y = (N + H)*cos(B)*sin(L);
pc.z = (N*(1 - e_*e_) + H)*sin(B);
return;
}

View File

@ -28,9 +28,8 @@
#ifndef _GCTL_REFELLIPSOID_H
#define _GCTL_REFELLIPSOID_H
#include "../core/macro.h"
#include "../core/exceptions.h"
#include "../maths/mathfunc.h"
#include "../poly/vertex.h"
#include "gmath.h"
namespace gctl
{
@ -38,6 +37,7 @@ namespace gctl
{
// reference system named as the planet name represents a mean radius reference sphere
Earth,
MagEarth,
Moon,
Mars,
// Reference ellipsoids
@ -75,17 +75,48 @@ namespace gctl
void set(double R, double r_or_flat, bool is_flat = false);
/**
* @brief Get the ellipsoidal radius at the inquiring latitude
* @brief Compute the local rRdius of curvature on the reference ellipsoid
*
* @param lati the inquiring latitude
* @return ellipsoidal radius
* @param geodetic_lati the inquiring latitude
* @return radius
*/
double radius_at(double lati);
double geodetic_radius(double geodetic_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_;
double r_, R_, f_; // semi-minor, semi-major and flat rate
double e_;
};
}
};
#endif // _GCTL_REFELLIPSOID_H

View File

@ -28,10 +28,7 @@
#ifndef _GCTL_SHAPEFUNC_H
#define _GCTL_SHAPEFUNC_H
#include "../gctl_config.h"
#include "../core/enum.h"
#include "../core/exceptions.h"
#include "../maths/mathfunc.h"
#include "gmath.h"
namespace gctl
{
@ -535,6 +532,6 @@ namespace gctl
size_t idx, eshape_value_e vt, edge_orient_e ot, double &xval, double &yval, double &zval,
double &xval2, double &yval2, double &zval2);
};
}
};
#endif // _GCTL_SHAPEFUNC_H

View File

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

Some files were not shown because too many files have changed in this diff Show More