tmp
This commit is contained in:
parent
11d29c3dd4
commit
4030726e33
@ -1,6 +1,6 @@
|
||||
cmake_minimum_required(VERSION 3.15.2)
|
||||
# 设置项目名称与语言
|
||||
project(GCTL VERSION 1.1)
|
||||
project(GCTL VERSION 2.0)
|
||||
# 添加配置配件编写的函数
|
||||
include(CMakePackageConfigHelpers)
|
||||
|
||||
|
||||
@ -132,7 +132,7 @@ namespace gctl
|
||||
|
||||
if (idx >= vec2d.size()) throw std::runtime_error("[gctl::str2type_vector2d] Invalid index.");
|
||||
|
||||
order_str = order_str.substr(order_str.find_first_of(',') + 1);
|
||||
order_str = order_str.substr(order_str.find_first_of(',') + 1);
|
||||
str2type_vector(vec2d[idx], vec);
|
||||
return;
|
||||
}
|
||||
|
||||
@ -28,8 +28,7 @@
|
||||
#ifndef _GCTL_MATH_H
|
||||
#define _GCTL_MATH_H
|
||||
|
||||
#include "math/mathfunc.h"
|
||||
#include "math/algorithm_func.h"
|
||||
#include "math/gmath.h"
|
||||
|
||||
#include "math/legendre.h"
|
||||
#include "math/linear_algebra.h"
|
||||
|
||||
@ -1,137 +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 "mathfunc.h"
|
||||
|
||||
namespace gctl
|
||||
{
|
||||
/**
|
||||
* @brief 按距离反比加权计算均值
|
||||
*
|
||||
* @param dis_vec 距离向量
|
||||
* @param val_vec 数值向量
|
||||
* @param[in] order 距离加权的阶次 默认为1
|
||||
*
|
||||
* @return 加权平均值
|
||||
*/
|
||||
double dist_inverse_weight(std::vector<double> *dis_vec, std::vector<double> *val_vec, int order = 1);
|
||||
|
||||
/**
|
||||
* @brief 查找一个数在已排序的数组中的位置，即找到包含该数的最小区间
|
||||
*
|
||||
* @param[in] in_array 输入数组
|
||||
* @param[in] array_size 数组大小
|
||||
* @param[in] in_val 查找值
|
||||
* @param index 返回的索引值
|
||||
*
|
||||
* @return 成功0失败-1
|
||||
*/
|
||||
int find_index(const double *in_array, int array_size, double in_val, int &index);
|
||||
|
||||
/**
|
||||
* @brief 查找一个数在已排序的数组中的位置，即找到包含该数的最小区间
|
||||
*
|
||||
* @param in_array 输入数组
|
||||
* @param[in] in_val 查找值
|
||||
* @param index 返回的索引值
|
||||
*
|
||||
* @return 成功0失败-1
|
||||
*/
|
||||
int find_index(array<double> *in_array, double in_val, int &index);
|
||||
|
||||
/**
|
||||
* @brief 计算一维分形模型
|
||||
*
|
||||
* @param out_arr 输出数组
|
||||
* @param[in] l_val 分形计算的左端点值
|
||||
* @param[in] r_val 分形计算的右端点值
|
||||
* @param[in] maxi_range 最大变化值
|
||||
* @param[in] smoothness 变化光滑度
|
||||
*/
|
||||
void fractal_model_1d(array<double> &out_arr, int out_size, double l_val,
|
||||
double r_val, double maxi_range, double smoothness);
|
||||
|
||||
/**
|
||||
* @brief 计算二维分形模型
|
||||
*
|
||||
* @param out_arr 输出数组
|
||||
* @param[in] dl_val 分形计算的左下角端点值
|
||||
* @param[in] dr_val 分形计算的右下角端点值
|
||||
* @param[in] ul_val 分形计算的左上角端点值
|
||||
* @param[in] ur_val 分形计算的右上角端点值
|
||||
* @param[in] maxi_range 最大变化值
|
||||
* @param[in] smoothness 变化光滑度
|
||||
*/
|
||||
void fractal_model_2d(_2d_matrix &out_arr, int r_size, int c_size, double dl_val,
|
||||
double dr_val, double ul_val, double ur_val, double maxi_range, double smoothness, unsigned int seed = 0);
|
||||
|
||||
/**
|
||||
* @brief 一维数组差分（使用二阶差分公式）
|
||||
*
|
||||
* 计算一个一维数组中相邻元素间的差分结果（求导）。
|
||||
*
|
||||
* @param[in] in 输入数组
|
||||
* @param diff 输出的差分结果
|
||||
* @param[in] spacing 相邻元素的距离
|
||||
* @param[in] order 求导的次数。最小为1（默认），最大为4。两边的数据将分别使用对应的向前或向后差分公式
|
||||
*/
|
||||
void difference_1d(const array<double> &in, array<double> &diff, double spacing, int order = 1);
|
||||
|
||||
/**
|
||||
* @brief 二维数组差分（使用二阶差分公式）
|
||||
*
|
||||
* 计算一个二维数组中相邻元素间的差分结果（求导）。
|
||||
*
|
||||
* @param[in] in 输入数组
|
||||
* @param diff 输出的差分结果
|
||||
* @param[in] spacing 相邻元素对应方向的距离
|
||||
* @param[in] d_type 求导的类型
|
||||
* @param[in] order 求导的次数。最小为1（默认），最大为4。边缘的数据将分别使用对应的向前或向后差分公式
|
||||
*/
|
||||
void difference_2d(const _2d_matrix &in, _2d_matrix &diff, double spacing, gradient_type_e d_type, int order = 1);
|
||||
|
||||
/**
|
||||
* @brief 二维数组差分（使用二阶差分公式）
|
||||
*
|
||||
* 计算一个二维数组中相邻元素间的差分结果（求导）。数组以列优先的方式储存为一个一维数组
|
||||
*
|
||||
* @param[in] in 输入数组
|
||||
* @param diff 输出的差分结果
|
||||
* @param[in] row_size 数组二维排列的行数
|
||||
* @param[in] col_size 数组二维排列的列数
|
||||
* @param[in] spacing 相邻元素对应方向的距离
|
||||
* @param[in] d_type 求导的类型
|
||||
* @param[in] order 求导的次数。最小为1（默认），最大为4。边缘的数据将分别使用对应的向前或向后差分公式
|
||||
*/
|
||||
void difference_2d(const array<double> &in, array<double> &diff, int row_size, int col_size,
|
||||
double spacing, gradient_type_e d_type, int order = 1);
|
||||
};
|
||||
|
||||
#endif // _GCTL_ALGORITHM_FUNC_H
|
||||
@ -28,7 +28,7 @@
|
||||
#ifndef _GCTL_EXTRAPOLATE_H
|
||||
#define _GCTL_EXTRAPOLATE_H
|
||||
|
||||
#include "mathfunc.h"
|
||||
#include "gmath.h"
|
||||
|
||||
namespace gctl
|
||||
{
|
||||
|
||||
@ -30,7 +30,7 @@
|
||||
|
||||
#include "../poly/triangle.h"
|
||||
#include "../poly/tetrahedron.h"
|
||||
#include "mathfunc.h"
|
||||
#include "gmath.h"
|
||||
|
||||
namespace gctl
|
||||
{
|
||||
|
||||
@ -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)
|
||||
{
|
||||
@ -25,15 +25,13 @@
|
||||
* 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
|
||||
|
||||
// system's head files
|
||||
#include <cmath>
|
||||
#include <random>
|
||||
#include <vector>
|
||||
|
||||
// library's head files
|
||||
#include "../core/macro.h"
|
||||
#include "../core/spmat.h"
|
||||
|
||||
@ -338,6 +336,107 @@ namespace gctl
|
||||
* @param[in] a_s 向量的个数
|
||||
*/
|
||||
void schmidt_orthogonal(const array<double> &a, array<double> &e, int a_s);
|
||||
|
||||
/**
|
||||
* @brief 按距离反比加权计算均值
|
||||
*
|
||||
* @param dis_vec 距离向量
|
||||
* @param val_vec 数值向量
|
||||
* @param[in] order 距离加权的阶次 默认为1
|
||||
*
|
||||
* @return 加权平均值
|
||||
*/
|
||||
double dist_inverse_weight(std::vector<double> *dis_vec, std::vector<double> *val_vec, int order = 1);
|
||||
|
||||
/**
|
||||
* @brief 查找一个数在已排序的数组中的位置，即找到包含该数的最小区间
|
||||
*
|
||||
* @param[in] in_array 输入数组
|
||||
* @param[in] array_size 数组大小
|
||||
* @param[in] in_val 查找值
|
||||
* @param index 返回的索引值
|
||||
*
|
||||
* @return 成功0失败-1
|
||||
*/
|
||||
int find_index(const double *in_array, int array_size, double in_val, int &index);
|
||||
|
||||
/**
|
||||
* @brief 查找一个数在已排序的数组中的位置，即找到包含该数的最小区间
|
||||
*
|
||||
* @param in_array 输入数组
|
||||
* @param[in] in_val 查找值
|
||||
* @param index 返回的索引值
|
||||
*
|
||||
* @return 成功0失败-1
|
||||
*/
|
||||
int find_index(array<double> *in_array, double in_val, int &index);
|
||||
|
||||
/**
|
||||
* @brief 计算一维分形模型
|
||||
*
|
||||
* @param out_arr 输出数组
|
||||
* @param[in] l_val 分形计算的左端点值
|
||||
* @param[in] r_val 分形计算的右端点值
|
||||
* @param[in] maxi_range 最大变化值
|
||||
* @param[in] smoothness 变化光滑度
|
||||
*/
|
||||
void fractal_model_1d(array<double> &out_arr, int out_size, double l_val,
|
||||
double r_val, double maxi_range, double smoothness);
|
||||
|
||||
/**
|
||||
* @brief 计算二维分形模型
|
||||
*
|
||||
* @param out_arr 输出数组
|
||||
* @param[in] dl_val 分形计算的左下角端点值
|
||||
* @param[in] dr_val 分形计算的右下角端点值
|
||||
* @param[in] ul_val 分形计算的左上角端点值
|
||||
* @param[in] ur_val 分形计算的右上角端点值
|
||||
* @param[in] maxi_range 最大变化值
|
||||
* @param[in] smoothness 变化光滑度
|
||||
*/
|
||||
void fractal_model_2d(_2d_matrix &out_arr, int r_size, int c_size, double dl_val,
|
||||
double dr_val, double ul_val, double ur_val, double maxi_range, double smoothness, unsigned int seed = 0);
|
||||
|
||||
/**
|
||||
* @brief 一维数组差分（使用二阶差分公式）
|
||||
*
|
||||
* 计算一个一维数组中相邻元素间的差分结果（求导）。
|
||||
*
|
||||
* @param[in] in 输入数组
|
||||
* @param diff 输出的差分结果
|
||||
* @param[in] spacing 相邻元素的距离
|
||||
* @param[in] order 求导的次数。最小为1（默认），最大为4。两边的数据将分别使用对应的向前或向后差分公式
|
||||
*/
|
||||
void difference_1d(const array<double> &in, array<double> &diff, double spacing, int order = 1);
|
||||
|
||||
/**
|
||||
* @brief 二维数组差分（使用二阶差分公式）
|
||||
*
|
||||
* 计算一个二维数组中相邻元素间的差分结果（求导）。
|
||||
*
|
||||
* @param[in] in 输入数组
|
||||
* @param diff 输出的差分结果
|
||||
* @param[in] spacing 相邻元素对应方向的距离
|
||||
* @param[in] d_type 求导的类型
|
||||
* @param[in] order 求导的次数。最小为1（默认），最大为4。边缘的数据将分别使用对应的向前或向后差分公式
|
||||
*/
|
||||
void difference_2d(const _2d_matrix &in, _2d_matrix &diff, double spacing, gradient_type_e d_type, int order = 1);
|
||||
|
||||
/**
|
||||
* @brief 二维数组差分（使用二阶差分公式）
|
||||
*
|
||||
* 计算一个二维数组中相邻元素间的差分结果（求导）。数组以列优先的方式储存为一个一维数组
|
||||
*
|
||||
* @param[in] in 输入数组
|
||||
* @param diff 输出的差分结果
|
||||
* @param[in] row_size 数组二维排列的行数
|
||||
* @param[in] col_size 数组二维排列的列数
|
||||
* @param[in] spacing 相邻元素对应方向的距离
|
||||
* @param[in] d_type 求导的类型
|
||||
* @param[in] order 求导的次数。最小为1（默认），最大为4。边缘的数据将分别使用对应的向前或向后差分公式
|
||||
*/
|
||||
void difference_2d(const array<double> &in, array<double> &diff, int row_size, int col_size,
|
||||
double spacing, gradient_type_e d_type, int order = 1);
|
||||
};
|
||||
|
||||
#endif // _GCTL_MATHFUNC_H
|
||||
#endif // _GCTL_GMATH_H
|
||||
@ -29,7 +29,6 @@
|
||||
#define _GCTL_INTERPOLATE_H
|
||||
|
||||
#include "geometry3d.h"
|
||||
#include "algorithm_func.h"
|
||||
|
||||
namespace gctl
|
||||
{
|
||||
|
||||
@ -1,504 +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 "mathfunc.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;
|
||||
}
|
||||
@ -29,7 +29,7 @@
|
||||
#define _GCTL_REFELLIPSOID_H
|
||||
|
||||
#include "../poly/vertex.h"
|
||||
#include "mathfunc.h"
|
||||
#include "gmath.h"
|
||||
|
||||
namespace gctl
|
||||
{
|
||||
|
||||
@ -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 "mathfunc.h"
|
||||
#include "gmath.h"
|
||||
|
||||
namespace gctl
|
||||
{
|
||||
|
||||
@ -30,6 +30,13 @@
|
||||
const int DaysInMonth[13] = {0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
|
||||
const std::string MonthName[13] = {"Null", "Jan.", "Feb.", "Mar.", "Apr."," May", "Jun.", "Jul.", "Aug.", "Sep.", "Oct.", "Nov.", "Dec."};
|
||||
|
||||
std::string gctl::time_stamp()
|
||||
{
|
||||
time_t now = time(0);
|
||||
char* dt = ctime(&now);
|
||||
return std::string(dt);
|
||||
}
|
||||
|
||||
// Function returns true if the given year is a leap year and returns false otherwise.
|
||||
bool gctl::is_leap_year(int year)
|
||||
{
|
||||
|
||||
@ -48,6 +48,9 @@ namespace gctl
|
||||
Year,
|
||||
};
|
||||
|
||||
// return a string of current time stamp.
|
||||
std::string time_stamp();
|
||||
|
||||
// Function returns true if the given year is a leap year and returns false otherwise.
|
||||
bool is_leap_year(int year);
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user