217 lines
7.8 KiB
C++
217 lines
7.8 KiB
C++
/********************************************************
|
||
* ██████╗ ██████╗████████╗██╗
|
||
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
||
* ██║ ███╗██║ ██║ ██║
|
||
* ██║ ██║██║ ██║ ██║
|
||
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
||
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
||
* 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_MATHFUNC_H
|
||
#define _GCTL_MATHFUNC_H
|
||
|
||
// 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
|
||
{
|
||
/**
|
||
* @brief Return a random number in the range [low, hig]
|
||
*
|
||
* @note Call srand(seed) to initiate the random squence before using this function.
|
||
*
|
||
* @param low Lower bound
|
||
* @param hig Higher bound
|
||
* @return Random value
|
||
*/
|
||
int random(int low, int hig);
|
||
|
||
/**
|
||
* @brief Return a random number in the range [low, hig]
|
||
*
|
||
* @note Call srand(seed) to initiate the random squence before using this function.
|
||
*
|
||
* @param low Lower bound
|
||
* @param hig Higher bound
|
||
* @return Random value
|
||
*/
|
||
double random(double low, double hig);
|
||
|
||
/**
|
||
* @brief Return a random number in the range [low, hig]
|
||
*
|
||
* @note Call srand(seed) to initiate the random squence before using this function.
|
||
*
|
||
* @param low Lower bound
|
||
* @param hig Higher bound
|
||
* @return Random value
|
||
*/
|
||
std::complex<double> random(std::complex<double> low, std::complex<double> hig);
|
||
|
||
/**
|
||
* @brief 比较两个浮点数。注意比较的精度为 eps 减一位，比如 eps 为1e-10则比较的精度为小数点后9位。
|
||
*
|
||
* @param f 第一个浮点数
|
||
* @param v 第二个浮点数
|
||
* @param eps 比较的精度
|
||
* @return true 两个浮点数相等
|
||
* @return false 两个浮点数不相等
|
||
*/
|
||
bool isequal(double f, double v, double eps = GCTL_ZERO);
|
||
|
||
/**
|
||
* @brief 符号函数
|
||
*
|
||
* @param[in] a 输入值
|
||
*
|
||
* @return 大于零返回1 否则返回-1 等于0则返回0
|
||
*/
|
||
double sign(double a);
|
||
|
||
/**
|
||
* @brief 二分法求一个正数的n次方根
|
||
*
|
||
* @param[in] val 输入值
|
||
* @param[in] order 开方次数
|
||
* @param[in] eps 终止精度
|
||
*
|
||
* @return 开方值
|
||
*/
|
||
double sqrtn(double val, int order, double eps = 1e-5);
|
||
|
||
/**
|
||
* @brief 计算经纬度网格的面积
|
||
*
|
||
* @param lon1 经度起点
|
||
* @param lon2 经度终点
|
||
* @param lat1 纬度起点
|
||
* @param lat2 纬度终点
|
||
* @param R 半径
|
||
* @return 面积
|
||
*/
|
||
double geographic_area(double lon1, double lon2, double lat1, double lat2, double R);
|
||
|
||
/**
|
||
* @brief 计算两个经纬度点之间的距离
|
||
*
|
||
* @param lon1 经度起点
|
||
* @param lon2 经度终点
|
||
* @param lat1 纬度起点
|
||
* @param lat2 纬度终点
|
||
* @param R 半径
|
||
* @return 面积
|
||
*/
|
||
double geographic_distance(double lon1, double lon2, double lat1, double lat2, double R);
|
||
|
||
/**
|
||
* @brief 计算一个椭圆在不同位置的半径
|
||
*
|
||
* @param[in] x_len x方向半径
|
||
* @param[in] y_len y方向半径
|
||
* @param[in] arc 计算方向与x轴正方向的夹角（弧度） 逆时针为正
|
||
* @param[in] x_arc x轴正方向绕原点逆时针旋转的角度 （弧度）
|
||
*
|
||
* @return 半径值
|
||
*/
|
||
double ellipse_radius_2d(double x_len, double y_len, double arc, double x_arc = 0.0);
|
||
|
||
/**
|
||
* @brief 计算已椭圆为基准面的高程数据的绝对坐标位置。
|
||
*
|
||
* 假设高程数据的测量方向（如大地水准面）为椭圆切线的垂直方向（与切点与坐标原点的连线方向不一致）
|
||
* 则高程点的绝对空间位置的维度值与切点的维度值也不一致。此函数可以计算校正后的维度值与球心半径。
|
||
*
|
||
* @param x_len 椭圆的x方向半径（一般为长轴）
|
||
* @param y_len 椭圆的y方向半径（一般为短轴）
|
||
* @param arc 计算方向与x轴正方向的夹角（弧度） 逆时针为正（等于纬度）
|
||
* @param elev 切点的高程值
|
||
* @param out_arc 校正后的维度值
|
||
* @param out_rad 校正后高程点的球心半径
|
||
* @param x_arc x轴正方向绕原点逆时针旋转的角度 （弧度），默认为0
|
||
*/
|
||
void ellipse_plus_elevation_2d(double x_len, double y_len, double arc, double elev,
|
||
double &out_arc, double &out_rad, double x_arc = 0.0);
|
||
|
||
/**
|
||
* @brief 椭球或者球在不同球面位置的半径
|
||
*
|
||
* @param[in] x_len x方向半径
|
||
* @param[in] y_len y方向半径
|
||
* @param[in] z_len z方向半径
|
||
* @param[in] phi 计算方向与x轴正方向的夹角（弧度） 逆时针为正
|
||
* @param[in] theta 计算方向与z轴正方向的夹角（弧度） 逆时针为正
|
||
*
|
||
* @return 半径值
|
||
*/
|
||
double ellipsoid_radius(double x_len, double y_len, double z_len, double phi, double theta);
|
||
|
||
/**
|
||
* @brief 使用牛顿迭代法计算一个矩阵的近似逆
|
||
*
|
||
* 此方法适合用于对角占优的方阵求逆。
|
||
*
|
||
* @param[in] in_mat 输入矩阵
|
||
* @param inverse_mat 逆矩阵
|
||
* @param[in] epsilon 迭代的终止精度
|
||
* @param[in] iter_times 迭代的次数。默认为10
|
||
* @param[in] initiate 对逆矩阵进行初始化
|
||
*
|
||
* @return 最终的迭代精度
|
||
*/
|
||
double newton_inverse(const _2d_matrix &in_mat, _2d_matrix &inverse_mat,
|
||
double epsilon = 1e-8, int iter_times = 10, bool initiate = true);
|
||
|
||
/**
|
||
* @brief 使用牛顿迭代法计算一个矩阵的近似逆
|
||
*
|
||
* 此方法适合用于对角占优的方阵求逆。
|
||
*
|
||
* @param[in] in_mat 输入矩阵（稀疏的）
|
||
* @param inverse_mat 逆矩阵
|
||
* @param[in] epsilon 迭代的终止精度
|
||
* @param[in] iter_times 迭代的次数。默认为10
|
||
* @param[in] initiate 对逆矩阵进行初始化
|
||
*
|
||
* @return 最终的迭代精度
|
||
*/
|
||
double newton_inverse(const spmat<double> &in_mat, _2d_matrix &inverse_mat,
|
||
double epsilon = 1e-8, int iter_times = 10, bool initiate = true);
|
||
|
||
/**
|
||
* @brief 使用施密特正交化将线性无关的向量组a转换为标准正交向量组e
|
||
*
|
||
* @note 经过正交变化的向量组中任意两个向量正交。
|
||
*
|
||
* @param[in] a 个数为a_s的线性无关的向量组的指针（这里将向量组保存为一个连续的一维数组）
|
||
* @param e 输出的转换后的标准正交向量组（保存为一个连续的一维数组）
|
||
* @param[in] a_s 向量的个数
|
||
*/
|
||
void schmidt_orthogonal(const array<double> &a, array<double> &e, int a_s);
|
||
}
|
||
|
||
#endif // _GCTL_MATHFUNC_H
|