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