sketchbook/common/math.hpp

#ifndef COMMON_MATH_HPP
#define COMMON_MATH_HPP
#include "simple/support.hpp"
#include "simple/geom.hpp"

namespace common
{

using namespace simple;

constexpr auto half = geom::vector<float,2>::one(.5f);
const float tau = 2*std::acos(-1);

// https://www.youtube.com/watch?v=tk58sBLWzHk
using proportion3 = geom::vector<float,3>;

proportion3 meet(proportion3 line1, proportion3 line2)
{
	const auto [x1, y1, z1] = line1;
	const auto [x2, y2, z2] = line2;
	// this be the mysterious cross product they say
	// the source above has the sing of z flipped, but that does not seem to
	// actually  matter as long as you're consistent
	return {y1*z2 - y2*z1, z1*x2 - z2*x1, x1*y2 - x2*y1};
}

proportion3 join(proportion3 point1, proportion3 point2)
{
	return meet(point1, point2);
}

proportion3 join(geom::vector<float,2> point1, geom::vector<float,2> point2)
{
	constexpr auto denominator = geom::vector(1.f);
	return join(point1.concat(denominator), point2.concat(denominator));
}

template <typename Vector>
[[nodiscard]] constexpr
Vector project_on_unit(Vector x, Vector surface)
{
	return surface * surface(x);
}

template <typename Vector>
[[nodiscard]] constexpr
Vector project(Vector x, Vector surface)
{
	return project_on_unit(x,surface) / surface.magnitude();
}

template <typename Vector>
[[nodiscard]] constexpr
Vector reject_from_unit(Vector x, Vector surface)
{
	return x - project_on_unit(x,surface);
}

template <typename Vector>
[[nodiscard]] constexpr
Vector reject(Vector x, Vector surface)
{
	return x - project(x,surface);
}

template <typename Vector>
[[nodiscard]] constexpr
Vector reflect_in_unit(Vector x, Vector surface)
{
	return 2*project_on_unit(x,surface) - x;
	// return project_on_unit(x,surface) - reject_from_unit(x,surface);
}

template <typename Vector>
[[nodiscard]] constexpr
Vector reflect(Vector x, Vector surface)
{
	return project(x,surface) - reject(x,surface);
}

template <typename Vector>
[[nodiscard]] constexpr
Vector rotate_scale(Vector x, Vector half_angle)
{
	x.y() = -x.y();
	return reflect_in_unit(x, half_angle);
	// return reflect_in_unit(
	// 		reflect_in_unit(x, Vector::i()),
	// 		half_angle );
}

template <typename Vector>
[[nodiscard]] constexpr
Vector rotate(Vector x, support::range<Vector> half_angle)
{
	return reflect( reflect(x, half_angle.lower()), half_angle.upper() );
}

template <typename Vector>
[[nodiscard]] constexpr
Vector rotate(Vector x, Vector half_angle)
{
	return rotate( x, support::range{Vector::i(), half_angle} );
}

template <typename Vector>
[[nodiscard]] constexpr
Vector normalize(Vector x)
{
	assert(x.quadrance() != 0);
	return x / support::root2(x.quadrance());
}


// see angle_from_cord_length.cpp
template <typename Number>
[[nodiscard]] constexpr
geom::vector<Number,2> cord_quadrance_angle(Number radius, Number cord_quadrance)
{
	auto mystery_ratio = cord_quadrance/(2*radius);
	return
	{
		radius - mystery_ratio,
		support::root2(cord_quadrance - mystery_ratio*mystery_ratio)
	};
}

template <typename Number>
[[nodiscard]] constexpr
geom::vector<Number,2> cord_length_angle(Number radius, Number cord_length)
{
	return cord_quadrance_angle(radius, cord_length * cord_length);
}

template <size_t Exponent = 3, typename Value = float>
class protractor
{
	public:
	using vector = geom::vector<Value,2>;

	using array = std::array<vector, (size_t(1)<<Exponent) + 1>;

	static constexpr array circle = []()
	{
		using support::halfway;

		array points{};

		points.front() = vector::i();
		points.back() = -vector::i();

		auto bisect = [](auto begin, auto end,
		auto& self)
		{
			if(end - begin <= 2)
				return;

			auto middle = halfway(begin, end);
			*(middle) = normalize(halfway(*begin, *(end - 1)));

			self(begin, middle + 1,
			self);
			self(middle, end,
			self);
		};

		if(Exponent > 0)
		{
			auto middle = halfway(points.begin(), points.end());
			*middle = vector::j();

			bisect(points.begin(), middle + 1,
			bisect);
			bisect(middle, points.end(),
			bisect);
		}

		return points;
	}();


	static constexpr vector tau(Value factor)
	{
		assert(Value{0} <= factor && factor <= Value{1});
		Value index = factor * (protractor::circle.size() - 1);
		int floor = index;
		Value fraction = index - floor;
		int ceil = floor + (fraction > 0);
		using support::way;
		return way(circle[floor], circle[ceil], fraction);
	}

	static constexpr Value tau(vector half_angle)
	{
		assert(half_angle.y() > 0);

		// too clever? I mean I could have done lower_bound upper_bound to be
		// oh so clever... this is just mildly clever... still fear that
		// find_if would optimize better
		auto lower = std::partition_point(circle.begin(), circle.end(),
			[half_angle](auto x)
			{
				auto side = rotate(x, tau(1/4.f));
				return side(half_angle) > 0;
			}
		);
		assert(lower != circle.end());
		auto upper = std::next(lower);

		auto index = lower - circle.begin();

		// NOTE: could just use half_angle itself if it was close enough to the circle, not sure if worth checking though
		auto on_circle = meet(join(vector::zero(),half_angle), join(lower,upper));

		auto sector_ratio = (on_circle - lower) / (upper - lower);
		return (index + sector_ratio.x()) / circle.size();
	}

	static constexpr vector tau()
	{
		return -vector::i();
	}

	static constexpr vector small_radian(Value value)
	{
		return {support::root2(Value{1} - value*value), value};
	}

	static vector radian(Value angle)
	{
		return {std::cos(angle), std::sin(angle)};
	}

	static constexpr vector angle(Value tau_factor)
	{
		if(tau_factor < (Value{1}/(circle.size()-1))/16 )
			return small_radian(tau_factor * common::tau);
		else
			return protractor::tau(tau_factor);
	}

};

float sinus(float angle)
{ return common::rotate(geom::vector<float,2>::i(), common::protractor<>::tau(support::wrap(angle,1.f))).y(); }

float cosinus(float angle)
{ return common::rotate(geom::vector<float,2>::i(), common::protractor<>::tau(support::wrap(angle,1.f))).x(); }

} // namespace common

#endif /* end of include guard */