#include "dl_math.h"

namespace dl::Math
{

std::vector<dl::PointU> GetCircle(const Size& size, const PointU& xy, unsigned rad)
{
	if (rad == 0)
		return { xy };

	std::vector<PointU> ret;

	PointU r = { rad, rad };
	PointU lt = Math::SubU(xy, r);
	PointU rb = Each::Limit(xy + r, {}, size);

	for (unsigned j = lt[1]; j <= rb[1]; ++j)
	{
		for (unsigned i = lt[0]; i <= rb[0]; ++i)
		{
			if (Math::DistancePow2U({i, j}, xy) <= rad * rad)
				ret.push_back(PointU{ i, j });
		}
	}
	return ret;
}

std::vector<dl::PointU> RingX::GetCircle(const Size& size, const PointU& xy, unsigned rad)
{
	if (rad == 0)
		return { xy };

	std::vector<PointU> ret;
	ret.reserve(GetCircleNumMax(rad));
	int r = (int)rad;

	Point xy_i = toPoint(xy);
	for (int y = -r; y <= r; ++y)
	{
		// 计算该行最大x偏移
		int maxX = (int)std::sqrt(r * r - y * y);
		for (int x = -maxX; x <= maxX; ++x) 
		{
			PointU p;
			if (Clamp(size, xy_i + Point{ x, y }, p))
				ret.push_back(p);
		}
	}
	assert(ret.size() <= GetCircleNumMax(rad));
	return ret;
}

std::vector<dl::PointU> RingX::Diffusion(const Size& size, const PointU& xy, unsigned count, FnDiffusionCheck func_valid, FnDiffusionCb func_cb)
{
	MapLimit map;
	map.SetSize(size);

	// 返回值
	std::vector<PointU> ret;
	// 是否已经遍历过
	std::unordered_map<unsigned, bool> map_index;

	// 当前需要扩散的格子
	std::list<PointU> list_open;
	// 加入首节点
	unsigned index = xy.x + xy.y * size.w;
	map_index[index] = true;
	list_open.push_back(xy);

	for (unsigned i = 0; i < count; ++i)
	{
		std::list<PointU> list_temp; // 此轮遍历到的格子
		for (auto& iter : list_open)
		{
			if (func_valid(iter))
			{
				// 添加到返回列表
				ret.push_back(iter);

				PointU xy_near;
				xy_near = iter;
				if (map.Left(xy_near))
				{
					unsigned index = xy_near.x + xy_near.y * size.w;
					if (!map_index.contains(index))
					{
						map_index[index] = true;
						list_temp.push_back(xy_near);
					}
				}

				xy_near = iter;
				if (map.Right(xy_near))
				{
					unsigned index = xy_near.x + xy_near.y * size.w;
					if (!map_index.contains(index))
					{
						map_index[index] = true;
						list_temp.push_back(xy_near);
					}
				}

				xy_near = iter;
				if (map.Bottom(xy_near))
				{
					unsigned index = xy_near.x + xy_near.y * size.w;
					if (!map_index.contains(index))
					{
						map_index[index] = true;
						list_temp.push_back(xy_near);
					}
				}

				xy_near = iter;
				if (map.Top(xy_near))
				{
					unsigned index = xy_near.x + xy_near.y * size.w;
					if (!map_index.contains(index))
					{
						map_index[index] = true;
						list_temp.push_back(xy_near);
					}
				}
			}
		}
		func_cb(i, list_temp);
		list_open = list_temp;
	}
	return ret;
}

std::vector<dl::Point> GetCircle(const Point& xy, unsigned rad)
{
	if (rad == 0)
		return { xy };

	std::vector<Point> ret;

	int r = rad;
	for (int j = -r; j <= r; ++j)
	{
		for (int i = -r; i <= r; ++i)
		{
			if (i * i + j * j <= r * r)
				ret.push_back(Point{ i, j } + xy);
		}
	}
	return ret;
}

dl::Position4 Slerp(const Position4& q0, const Position4& q1, const Float& t)
{
	Float magnitude = sqrt(DistancePow2(q0) * DistancePow2(q1));
	assert(magnitude > 0);
	Float product = DotProduct(q0, q1) / magnitude;
	Float product_abs = std::abs(product);
	if (product_abs < 1_f - std::numeric_limits<Float>::epsilon())
	{
		Float theta = acos(product_abs);
		Float d = sin(theta);
		assert(d > 0);

		Float sign = product < 0 ? -1_f : 1_f;
		Float s0 = sin((1_f - t) * theta) / d;
		Float s1 = sin(sign * t * theta) / d;

		return q0 * s0 + q1 * s1;
	}
	else
		return q0;
}

std::vector<dl::Point> DrawLineLerp(const Point& p0, const Point& p1)
{
	std::vector<Point> ret;
	int N = DistanceDiagonal(p0, p1);
	for (int step = 0; step <= N; ++step)
	{
		double t = N ? ((double)step / N) : 0.0;
		ret.push_back(Each::Lerp(p0, p1, t, std::round));
	}
	return ret;
}

std::vector<dl::Point> DrawLineBresenham(const Point& p0, const Point& p1)
{
	std::vector<Point> ret;

	const int& x0 = p0[0];
	const int& y0 = p0[1];
	const int& x1 = p1[0];
	const int& y1 = p1[1];

	int dx = x1 - x0;
	int dy = y1 - y0;
	int ux = dx > 0 ? 1 : -1;
	int uy = dy > 0 ? 1 : -1;
	dx = abs(dx);
	dy = abs(dy);
	int dx_2 = dx * 2;
	int dy_2 = dy * 2;
	if (dx > dy)
	{
		int e = -dx;
		int x = x0;
		int y = y0;

		for (x = x0; x != x1; x += ux)
		{
			ret.push_back({ x, y });
			e = e + dy_2;
			if (e > 0)
			{
				y += uy;
				e = e - dx_2;
			}
		}
		ret.push_back({ x, y });
	}
	else
	{
		int e = -dy;
		int x = x0;
		int y = y0;
		for (y = y0; y != y1; y += uy)
		{
			ret.push_back({ x, y });
			e = e + dx_2;
			if (e > 0)
			{
				x += ux;
				e = e - dy_2;
			}
		}
		ret.push_back({ x, y });
	}
	return ret;
}

std::vector<dl::Point> DrawLineOrth(const Point& p0, const Point& p1)
{
	std::vector<Point> ret;

	int dx = p1[0] - p0[0];
	int dy = p1[1] - p0[1];
	int nx = abs(dx);
	int ny = abs(dy);
	int sign_x = dx > 0 ? 1 : -1;
	int sign_y = dy > 0 ? 1 : -1;

	Point p = p0;
	ret.push_back(p);
	int ix = 0;
	int iy = 0;
	for (; ix < nx || iy < ny;)
	{
		int decision = (1 + 2 * ix) * ny - (1 + 2 * iy) * nx;
		if (decision < 0)
		{//水平
			p[0] += sign_x;
			ix++;
		}
		else
		{//垂直
			p[1] += sign_y;
			iy++;
		}
		ret.push_back(p);
	}
	return ret;
}

std::vector<dl::Point> DrawLineCover(const Point& p0, const Point& p1)
{
	std::vector<Point> ret;

	int dx = p1[0] - p0[0];
	int dy = p1[1] - p0[1];
	int nx = abs(dx);
	int ny = abs(dy);
	int sign_x = dx > 0 ? 1 : -1;
	int sign_y = dy > 0 ? 1 : -1;

	Point p = p0;
	ret.push_back(p);
	int ix = 0;
	int iy = 0;
	for (; ix < nx || iy < ny;)
	{
		int decision = (1 + 2 * ix) * ny - (1 + 2 * iy) * nx;
		if (decision == 0)
		{//对角线
			p[0] += sign_x;
			p[1] += sign_y;
			ix++;
			iy++;
		}
		else if (decision < 0)
		{//水平
			p[0] += sign_x;
			ix++;
		}
		else
		{//垂直
			p[1] += sign_y;
			iy++;
		}
		ret.push_back(p);
	}
	return ret;
}

void MergePath(std::list<Point>& list_pos)
{
	size_t size = list_pos.size();
	if (list_pos.size() < 3)
		return;//少于3个节点不用处理
	//初始向量
	auto iter_0 = list_pos.begin();//前一个节点
	auto iter_1 = iter_0;//后一个节点
	++iter_1;
	Point v = *iter_1 - *iter_0;//向量指向前一个节点
	++iter_0;
	++iter_1;

	//直到后一个到end
	for (; iter_1 != list_pos.end();)
	{
		Point v1 = *iter_1 - *iter_0;
		if (v1 == v)
		{//方向相同，删除前一个节点
			list_pos.erase(iter_0);
			iter_0 = iter_1++;//iter_0失效、指向iter_1，iter_1指向下一位置
		}
		else
		{
			v = v1;
			++iter_0;
			++iter_1;
		}
	}
}

double GetSolarDeclination(int day)
{
	// 近似公式：赤纬 = 23.44° × sin(360/365 × (284 + N))
	// 这是地理学中常用的经验公式
	double angle = (360.0 / 365.0) * (284 + day);
	double declination = 23.44 * std::sin(DegreesToRadians(angle));
	return declination;
}

double GetSolarElevation(double latitude, double declination)
{
	// H = 90° - |纬度 - 赤纬|
	double diff = fabs(latitude - declination);
	double elevation = 90.0 - diff;
	return elevation;
}

int GetSeason(int day, double latitude)
{
	bool is_northern = (latitude >= 0);

	// 定义节气的大致日期（北半球）
	// 春分：3月20日（第79天）
	// 夏至：6月21日（第172天）
	// 秋分：9月23日（第266天）
	// 冬至：12月22日（第356天）

	int spring_equinox = 79;
	int summer_solstice = 172;
	int autumn_equinox = 266;
	int winter_solstice = 356;

	if (!is_northern)
	{
		// 南半球：季节相反，交换春分/秋分，夏至/冬至
		std::swap(spring_equinox, autumn_equinox);
		std::swap(summer_solstice, winter_solstice);
	}

	// 判断季节（使用环形日期）
	// 为了处理跨年，将日期映射到0-364的环上
	int day_ring = day % 365;

	// 简化判断：将四个节点之间的区间分配给四个季节
	// 这里用区间中点划分更准确
	if (day_ring >= spring_equinox && day_ring < summer_solstice)
	{
		return 0;  // 春季
	}
	else if (day_ring >= summer_solstice && day_ring < autumn_equinox)
	{
		return 1;  // 夏季
	}
	else if (day_ring >= autumn_equinox && day_ring < winter_solstice)
	{
		return 2;  // 秋季
	}
	else
	{
		return 3;  // 冬季（包含跨年部分）
	}
}

double GetBaseTemperature(double latitude, int day)
{
	// 基础纬度降温效应（年均温基准）
	double base_temp = 26.0 - (std::fabs(latitude) * 0.6);

	// 太阳高度角
	double solar_elevation = GetSolarElevation(latitude, GetSolarDeclination(day));

	// 太阳高度角修正系数
	// 高度角为90°时热力系数为1.0，为0°时系数为0.0
	double heat_factor = std::sin(DegreesToRadians(solar_elevation));

	// 确保热力系数在合理范围
	heat_factor = (std::max)(0.1, (std::min)(1.0, heat_factor));

	// 季节修正（太阳高度角已经反映了季节变化，但这里做额外微调）
	double seasonal_offset = 0.0;
	if (GetSeason(day, latitude) == 1)
	{
		seasonal_offset = 2.0;  // 夏季额外增温
	}
	else if (GetSeason(day, latitude) == 3)
	{
		seasonal_offset = -2.0; // 冬季额外降温
	}

	// 4. 组合计算：基准温度 × 热力系数 + 季节偏移
	double estimated_temp = base_temp * heat_factor + seasonal_offset;

	return estimated_temp;
}

}