so now i am trying to port it to c++ to organize it better this is harder, i can't do c++ effectively while confused on an ipad i made the attached c++ code today at a student center that i think or hope is also open to the public. it performs the same code for the inverse involute, calculates optimized points for an involute curve rasterized to 1.0 unit output precision via subdivision, and starts to lay out classes for a gear assembly. it's another stretch for me to use c++ again. excited the rasterization seems to work.
#include <numbers> #include <cmath> #include <cstdlib>
#include <vector>
#include <iostream>
typedef long double scalar;
typedef long long integer;
constexpr auto pi = std::numbers::pi_v<scalar>;
inline scalar powint(scalar base, unsigned int const & exponent)
{ return exponent > 0 ? base * powint(base, exponent-1) : 1; }
struct Rational
{
// i think i understand better if the denominator is kept smaller
integer numerator, denominator;
template <typename T> operator T const() const
{ return std::div(numerator, denominator); }
scalar to_scalar() const;
scalar operator/(scalar const & other);
integer max_units() const
{ return 1 + (numerator - 1) / denominator; }
integer max_per_unit() const
{ return 1 + (denominator - 1) / numerator; }
};
struct Point {
scalar x, y;
scalar distance(Point & other) const
{ return std::sqrt(powint(x-other.x,2) + powint(y-other.y,2)); }
bool within_cartesian_precision_of(Rational precision, Point const & other) const;
};
struct GearAssembly
{
scalar const angle_pressure;
scalar const radius_pitch_min;
int const tooth_count_min;
scalar circumference_pitch_min() const
{ return 2*pi*radius_pitch_min; }
scalar tooth_arcwidth_pitch() const
{ return circumference_pitch_min()/tooth_count_min; }
};
struct GearProfile
{
GearAssembly const & assembly;
int const tooth_count;
Point center;
scalar angle;
scalar circumference_pitch() const
{ return tooth_count * assembly.tooth_arcwidth_pitch(); }
scalar radius_pitch() const
{ return circumference_pitch()/(2*pi); }
scalar involute_radius_base() const
{ return radius_pitch() * std::cos(assembly.angle_pressure); }
void draw_gear() {
// draft
}
};
struct Involute
{
Point center_base;
scalar radius_base;
scalar angle_base;
static scalar pressure_to_involute(scalar angle)
{ return std::tan(angle) - angle; }
static long double involute_to_pressure(long double angle);
static void test_involute_to_pressure();
static scalar involute_to_rolling(scalar angle)
{ return angle + involute_to_pressure(angle); }
static scalar rolling_to_involute(scalar angle)
{ return angle - std::atan(angle); }
Point operator()(scalar angle_rolling) const
{
return {
center_base.x + radius_base * (
std::cos(angle_base + angle_rolling)
+ angle_rolling * std::sin(angle_base + angle_rolling)
),
center_base.y + radius_base * (
std::sin(angle_base + angle_rolling)
- angle_rolling * std::cos(angle_base + angle_rolling)
)
};
}
std::vector<Point> & points(
scalar angle_involute_start, scalar angle_involute_end,
Rational precision, std::vector<Point> & output
)
{
scalar dir = angle_involute_start < angle_involute_end ? 1 : -1;
struct next_data {
scalar angle;
Point point;
};
std::vector <next_data> later_points;
scalar angle_rolling_start = involute_to_rolling(angle_involute_start);
scalar angle_rolling_end = involute_to_rolling(angle_involute_end);
later_points.emplace_back(angle_rolling_end, (*this)(involute_to_rolling(angle_involute_end)));
scalar current_rolling = angle_rolling_start;
output.emplace_back((*this)(current_rolling));
Point * current_point = & output.back();
scalar next_rolling, length, projected_length;
Point next_point, line_vector, projected_point;
while (later_points.size()) {
next_rolling = (current_rolling + later_points.back().angle) / 2;
next_point = (*this)(next_rolling);
line_vector = later_points.back().point; line_vector.x -= current_point->x; line_vector.y -= current_point->y;
scalar length = std::sqrt(powint(line_vector.x,2) + powint(line_vector.y,2));
line_vector.x /= length; line_vector.y /= length;
projected_length =
line_vector.x * (next_point.x - current_point->x) +
line_vector.y * (next_point.y - current_point->y);
projected_point.x = line_vector.x * projected_length + current_point->x;
projected_point.y = line_vector.y * projected_length + current_point->y;
if (!next_point.within_cartesian_precision_of(precision, projected_point)) {
later_points.emplace_back(next_rolling, next_point);
continue;
} else {
output.emplace_back(std::move(later_points.back().point));
current_point = & output.back();
current_rolling = later_points.back().angle;
later_points.pop_back();
}
}
return output;
}
};
#include <iostream>
int main()
{
Involute::test_involute_to_pressure();
Involute example_curve_from_X{{10,10},{10},0};
std::vector<Point> points;
for (auto & point : example_curve_from_X.points(0, 90*pi/180, {1, 1}, points))
{
std::cout << "(" << point.x << "," << point.y << "),";
}
std::cout << std::endl;
}
// i'm not certain of the approaches
// concerns include overflow converting integer to scalar, as well as inaccuracy using scalar in divisions
// the goals below may be to perform as many of the operations using integer-with-float as possible, to reduce overflow,
// as well as to delay division
scalar Rational::to_scalar() const {
return std::abs(numerator) > std::abs(denominator) ?
numerator / scalar(denominator) :
scalar(numerator) / denominator;
}
scalar Rational::operator/(scalar const & other)
{
return std::abs(numerator) > std::abs(denominator) ?
numerator / (other * denominator) :
(numerator / other) / denominator;
}
bool Point::within_cartesian_precision_of(Rational precision, Point const & other) const
{
return int((this->x-0.5) * precision.denominator)/precision.numerator == int((other.x-0.5) * precision.denominator)/precision.numerator &&
int((this->y-0.5) * precision.denominator)/precision.numerator == int((other.y-0.5) * precision.denominator)/precision.numerator;
}
long double Involute::involute_to_pressure(long double angle) {
constexpr auto pi = std::numbers::pi_v<decltype(angle)>;
decltype(angle) inv;
if (-1.8 <= angle && angle <= 1.8) {
inv =
std::cbrt(3*angle)
- 2*angle/5
+ 9*std::cbrt(9*powint(angle,5))/175
- 2*std::cbrt(3*powint(angle,7))/175
- 144*std::cbrt(powint(angle,9))/67375
+ 3258*std::cbrt(9*powint(angle,11))/3128125
- 49711*std::cbrt(3*powint(angle,13))/153278125
- 1130112*std::cbrt(powint(angle,15))/9306171875
//+ 5169659643*std::cbrt(9*powint(angle,17))/95304506171875
//...
;
} else if (-5 < angle && angle < 5) {
auto zeta = angle - (std::tan(5*pi/12) - 5*pi/12);
inv =
5*pi/12
+ (7 - 4*std::sqrt(3))*zeta
- (388 - 224*std::sqrt(3))*powint(zeta,2)
+ (41784 - 24124*std::sqrt(3))*powint(zeta,3)
- (5623484 - 3246720*std::sqrt(3))*powint(zeta,4)
+ (847344640 - 489214656*std::sqrt(3))*powint(zeta,5)
- (410335269916.L/3 - 236907178544*std::sqrt(3)/3)*powint(zeta,6)
//+ (69385871593784.L/3 - 40059951642628*std::sqrt(3)/3)*powint(zeta,7)
//- (12132384030575012.L/3 - 7004635185964400*std::sqrt(3)/3)*powint(zeta,8)
;
} else {
inv =
pi/2
- 1/angle
+ pi/(2*powint(angle,2))
- (powint(pi/2,2) + 2/3)/powint(angle,3)
+ (powint(pi/2,3) + pi)/powint(angle,4)
- (powint(pi/2,4) + powint(pi,2) + 13/15)/powint(angle,5)
// + ...
;
}
//for (int i = 0; i < 4; ++ i) {
// auto tan_inv = std::tan(inv);
// inv += (angle - (tan_inv - inv))/(tan_inv * tan_inv);
//}
return inv;
}
#include <iostream>
void Involute::test_involute_to_pressure()
{
for (scalar u_deg = 5; u_deg < 90; u_deg += 5) {
scalar u = u_deg * pi / 180;
scalar x = pressure_to_involute(u);
scalar u2 = involute_to_pressure(x);
std::cout << u_deg << " " << x << " " << std::abs(u-u2) << std::endl;
}
}
makefile
Description: Binary data
