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

#include "point.h"

namespace gctl
{
	typedef point3s<double> point3ds;
	typedef point3s<float>  point3fs;

	/**
	 * @brief      球坐标系中的一个点（经纬度表示的）
	 */
	template <typename T>
	struct point3s
	{
		T rad, lon, lat;
		// 构造函数与拷贝构造函数
		point3s();
		point3s(T r, T ln, T lt);
		point3s(const point3s<T> &b);
		// 析构函数
		virtual ~point3s(){}
		// 结构体方法
		bool valid() const;
		void set(T r, T ln, T lt);
		void set(const point3s<T> &b);
		double module() const;
		point3c<T> s2c() const; // 转换为球坐标点
		void str(std::string str);
		/**
		 * @brief      Sets the i/o precision of the type.
		 *
		 * @param[in]  psn   The desired precision
		 */
		void set_io_precision(int psn);
		/**
		 * @brief      输出位置
		 *
		 * @param      os    输出流
		 */
		void out_loc(std::ostream &os, char deli) const;
		/**
		 * @brief      输入位置
		 *
		 * @param      os    输入流
		 */
		void in_loc(std::istream &os);
		/**
		 * @brief      返回位置
		 *
		 * @return     位置
		 */
		point3s<T> get_loc() const;
		/**
		 * @brief 返回正常球坐标中的Phi值
		 * 
		 * @return double 
		 */
		T get_phi() const;
		/**
		 * @brief 返回正常球坐标中的Theta值
		 * 
		 * @return double 
		 */
		T get_theta() const;
	};

	template <typename T>
	gctl::point3s<T>::point3s()
	{
		rad = lon = lat = NAN;
	}

	template <typename T>
	gctl::point3s<T>::point3s(T r, T ln, T lt)
	{
		set(r, ln, lt);
	}

	template <typename T>
	gctl::point3s<T>::point3s(const point3s<T> &b)
	{
		set(b);
	}

	template <typename T>
	bool gctl::point3s<T>::valid() const
	{
		if (std::isnan(rad) || std::isnan(lon) || std::isnan(lat)) return false;
		if (std::isinf(rad) || std::isinf(lon) || std::isinf(lat)) return false;
		if (rad < 0 || lon > 180 || lon < -180 || lat > 90 || lat < -90) return false;
		return true;
	}

	template <typename T>
	void gctl::point3s<T>::set(T r, T ln, T lt)
	{
		if (std::isnan(r) || std::isnan(ln) || std::isnan(lt) || 
			std::isinf(r) || std::isinf(ln) || std::isinf(lt))
		{
			throw invalid_argument("Invalid value detected. From point3ds::set(...)");
		}

		if (r < 0 || ln > 180 || ln < -180 || lt > 90 || lt < -90)
		{
			throw out_of_range("Invalid values. From point3s::set(...)");
		}

		rad = r; lon = ln; lat = lt;
		return;
	}

	template <typename T>
	void gctl::point3s<T>::set(const point3s<T> &b)
	{
		if (std::isnan(b.rad) || std::isnan(b.lon) || std::isnan(b.lat) || 
			std::isinf(b.rad) || std::isinf(b.lon) || std::isinf(b.lat))
		{
			throw invalid_argument("Invalid value detected. From point3s::set(...)");
		}

		if (b.rad < 0 || b.lon > 180 || b.lon < -180 || b.lat > 90 || b.lat < -90)
		{
			throw out_of_range("Invalid radius. From point3s::set(...)");
		}

		rad = b.rad; lon = b.lon; lat = b.lat;
		return;
	}

	template <typename T>
	double gctl::point3s<T>::module() const
	{
		return rad;
	}

	template <typename T>
	gctl::point3c<T> gctl::point3s<T>::s2c() const
	{
		/*
		point3c<T> outc;
		outc.x = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*cos((2.0 + lon/180.0)*GCTL_Pi);
		outc.y = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*sin((2.0 + lon/180.0)*GCTL_Pi);
		outc.z = rad*cos((0.5 - lat/180.0)*GCTL_Pi);
		return outc;
		*/

		point3c<T> v;
		v.x = rad*cos(GCTL_Pi*lat/180.0)*cos(GCTL_Pi*lon/180.0);
		v.y = rad*cos(GCTL_Pi*lat/180.0)*sin(GCTL_Pi*lon/180.0);
		v.z = rad*sin(GCTL_Pi*lat/180.0);
		return v;
	}

	template <typename T>
	void gctl::point3s<T>::str(std::string str)
	{
		std::smatch ret;
		std::regex pattern("\\((-?\\d*\\.?\\d+?),[ ]*(-?\\d*\\.?\\d+?),[ ]*(-?\\d*\\.?\\d+?)\\)");

		if (regex_search(str, ret, pattern))
		{
			rad = atof(std::string(ret[1]).c_str());
			lon = atof(std::string(ret[2]).c_str());
			lat = atof(std::string(ret[3]).c_str());
			return;
		}

		throw runtime_error("Fail to parse the input string: " + str + ". From point3s::str(...)");
	}

	template <typename T>
	void gctl::point3s<T>::set_io_precision(int psn)
	{
		if (psn < 0)
		{
			throw invalid_argument("Invalid precision. From point3s::set_io_precision(...)");
		}

		io_psn = psn;
		return;
	}

	template <typename T>
	void gctl::point3s<T>::out_loc(std::ostream &os, char deli) const
	{
		os << std::setprecision(io_psn) << rad << deli << lon << deli << lat;
		return;
	}

	template <typename T>
	void gctl::point3s<T>::in_loc(std::istream &os)
	{
		os >> rad >> lon >> lat;
		return;
	}

	template <typename T>
	gctl::point3s<T> gctl::point3s<T>::get_loc() const
	{
		return point3s<T>(rad, lon, lat);
	}

	template <typename T>
	T gctl::point3s<T>::get_phi() const
	{
		if (lon < .0) return GCTL_Pi*(2.0 - lon/180.0);
		return GCTL_Pi*lon/180.0; 
	}

	template <typename T>
	T gctl::point3s<T>::get_theta() const
	{
		return GCTL_Pi*(0.5 - lat/180.0);
	}

	template <typename T>
	bool operator ==(const point3s<T> &a, const point3s<T> &b)
	{
		if( fabs(a.lon - b.lon) < GCTL_ZERO && 
			fabs(a.lat - b.lat) < GCTL_ZERO && 
			fabs(a.rad - b.rad) < GCTL_ZERO)
		{
			return true;
		}
		else return false;
	}

	template <typename T>
	bool operator !=(const point3s<T> &a, const point3s<T> &b)
	{
		if( fabs(a.lon - b.lon) >= GCTL_ZERO || 
			fabs(a.lat - b.lat) >= GCTL_ZERO || 
			fabs(a.rad - b.rad) >= GCTL_ZERO)
		{
			return true;
		}
		else return false;
	}

	template <typename T>
	point3s<T> operator -(const point3s<T> &a, const point3s<T> &b)
	{
		point3s<T> m;
		m.rad = a.rad - b.rad;
		m.lon = a.lon - b.lon;
		m.lat = a.lat - b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator +(const point3s<T> &a, const point3s<T> &b)
	{
		gctl::point3s<T> m;
		m.rad = a.rad + b.rad;
		m.lon = a.lon + b.lon;
		m.lat = a.lat + b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(int sign, const point3s<T> &b)
	{
		gctl::point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(float sign, const point3s<T> &b)
	{
		gctl::point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(double sign, const point3s<T> &b)
	{
		gctl::point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(const point3s<T> &b, int sign)
	{
		gctl::point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(const point3s<T> &b, float sign)
	{
		gctl::point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(const point3s<T> &b, double sign)
	{
		point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	std::ostream &operator <<(std::ostream & os, const point3s<T> &a)
	{
		os << std::setprecision(io_psn) << a.rad << " " << a.lon << " " << a.lat;
		return os;
	}

	template <typename T>
	std::istream &operator >>(std::istream & os, point3s<T> &a)
	{
		os >> a.rad >> a.lon >> a.lat;
		return os;
	}

	/**
	 * @brief      函数判断两个point3ds类型是否等于
	 *
	 * @param[in]  a     三维空间内的一个向量或实点a。
	 * @param[in]  b     三维空间内的一个向量或实点b。
	 * @param[in]  cut_off  截断误差，默认值为 gctl_macro.h 中预设的 GCTL_ZERO 值，可在调用时改变。
	 * 注意比较的精度为 cut_off 减一位，比如 cut_off 为1e-10则比较的精度为小数点后9位。
	 *
	 * @return     是否相等。
	 */
	template <typename T>
	bool isequal(const point3s<T> &a, const point3s<T> &b, double cut_off = GCTL_ZERO)
	{
		if( fabs(a.rad - b.rad) < cut_off && 
			fabs(a.lon - b.lon) < cut_off && 
			fabs(a.lat - b.lat) < cut_off)
		{
			return true;
		}
		else return false;
	}

	/**
	 * @brief      函数判断两个point3ds类型是否不等于
	 *
	 * @param[in]  a     三维空间内的一个向量或实点a。
	 * @param[in]  b     三维空间内的一个向量或实点b。
	 * @param[in]  cut_off  截断误差，默认值为 gctl_macro.h 中预设的 GCTL_ZERO 值，可在调用时改变。
	 *
	 * @return     是否不相等。
	 */
	template <typename T>
	bool notequal(const point3s<T> &a, const point3s<T> &b, double cut_off = GCTL_ZERO)
	{
		if( fabs(a.rad - b.rad) >= cut_off || 
			fabs(a.lon - b.lon) >= cut_off || 
			fabs(a.lat - b.lat) >= cut_off)
		{
			return true;
		}
		else return false;
	}

	/**
	 * @brief      生成直球坐标系的点数组
	 *
	 * @param      out_ps  返回的点数组
	 * @param[in]  lonmin  lon最小值
	 * @param[in]  lonmax  lon最大值
	 * @param[in]  latmin  lat最小值
	 * @param[in]  latmax  lat最大值
	 * @param[in]  dlon    lon间隔
	 * @param[in]  dlat    lat间隔
	 * @param[in]  rad     半径值
	 * @param[in]  cor     网格起始位置
	 */
	template <typename T>
	void grid_points_2d(array<point3s<T>> &out_ps, T lonmin, T lonmax, T latmin, 
		T latmax, T dlon, T dlat, T rad, matrix_corner_e cor = BtmLeft)
	{
		if (lonmin >= lonmax || latmin >= latmax || 
			lonmin + dlon > lonmax || latmin + dlat > latmax || 
			dlon <= 0 || dlat <= 0)
		{
			throw invalid_argument("[gctl::grid_points_2d] Invalid parameters.");
		}

		int lonnum = round((lonmax-lonmin)/dlon) + 1;
		int latnum = round((latmax-latmin)/dlat) + 1;

		out_ps.resize(lonnum*latnum);
		if (cor == BtmLeft)
		{
			for (int j = 0; j < latnum; j++)
			{
				for (int i = 0; i < lonnum; i++)
				{
					out_ps.at(j*lonnum+i).lon = lonmin + dlon*i;
					out_ps.at(j*lonnum+i).lat = latmin + dlat*j;
					out_ps.at(j*lonnum+i).rad = rad;
				}
			}
		}
		else // cor == TopLeft
		{
			for (int j = 0; j < latnum; j++)
			{
				for (int i = 0; i < lonnum; i++)
				{
					out_ps.at(j*lonnum+i).lon = lonmin + dlon*i;
					out_ps.at(j*lonnum+i).lat = latmax - dlat*j;
					out_ps.at(j*lonnum+i).rad = rad;
				}
			}
		}
		return;
	}
};

#endif // _GCTL_POINT3S_H