gctl/lib/math/mathfunc.h

/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * 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

// system's head files
#include <cmath>
#include <random>
#include <vector>

// library's head files
#include "../core/macro.h"
#include "../core/spmat.h"

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]
     * 
     * @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