[BUG] Gears with high pressure angle self-intersect

[xloem]

Describe the bug
Gears with high pressure angles produce self-intersecting geometry. The implementation code appears to assume the pressure angle is low enough that the involute curves do not intersect each other before reaching the root radius.

An example of gears with these high pressure angles is the popular bearing at https://www.thingiverse.com/thing:53451 from 2013 by emmett.

Code To Reproduce Bug

include <BOSL2/std.scad>
include <BOSL2/gears.scad>

spur_gear(teeth=4,pressure_angle=60);

Expected behavior
Self intersections can produce problems with rendering and use. It is much better to have a final shape that is clearly closed.

Screenshots

Additional context
Add any other context about the problem here.

• OpenSCAD Version: This is a build of git commit 3241bf88b from 2023-01-04 . I'll be pulling a new version today, but it will take me a while to build and use it.
• Other libraries used: BOSL2 0f09d1f .

This problem has been solved a handful of ways elsewhere. I handwrote some gear code in C++ once that used the inverse of the involute function to calculate the intersection point here. Because this inverse is hard to calculate I attach this code:
involute_inverse.cpp . I am still new to OpenSCAD.

[adrianVmariano]

Note, the OpenSCAD version won't make a difference. (You can download prebuilt nightly builds unless you have special needs.)

We are absolutely NOT experts on gears, so insight you have on how to improve the gear implementation is welcome. One of our main open issues is that worm gears self-intersect. I didn't even know anybody ever used pressure angles above 20 degrees.

The OpenSCAD language is quirky, to say the least, but we can likely assist you if you want to make specific changes, for example. Revar wrote the core of the gear code and he's not very available these days. It looks like it's got a bunch of linear interpolation inverse lookups. Perhaps it can/should be completely rewritten? (This is the function _gear_tooth_profile() I'm referring to.)

Testing suggests that pressure angles above about 33 fail. Does it suffice to terminate the involute curve at a higher point, or are there other complications at the high pressure angles?

That inverse involute seems rather complicated indeed. It looks like a many term series approximation to the inverse function. I suspect that if we only need a single inverse we may find it adequate to use root_find(), our generic root finder, but I suppose the explicit series could be more robust.

[adrianVmariano]

f=function (angle) tan(angle) - angle*PI/180;

function inv_f(val) = root_find(function(x) f(x)-val, 0, 89.99);

for(u=[5:5:89])
{
  fu = f(u);
  rt=inv_f(fu);
  echo(u, fu, rt, rt-u);
}

Above runs the same test as your C++ code. Error is higher mostly because we don't have long doubles, I think. Translation of the other code to OpenSCAD if I didn't botch it:

  function involute_to_pressure(angle) = 
    angle==0 ? 0
  :  let(
        inv = abs(angle)<=1.8 ?
                  let(
                        angle1_3 = angle^(1/3),
                        angvec = [angle1_3, angle, angle1_3^5, angle1_3^7, angle1_3^9, angle1_3^11, angle1_3^13, angle1_3^15],
                        coef = [3^(1/3), -2/5, 9*9^(1/3)/175, -2*3^(1/3)/175, -144/67375, 3258*9^(1/3)/3128125, -49711*3^(1/3)/153278125, -1130112/9306171875]
                  )
                  coef*angvec
            : abs(angle)<5 ?
                  let(
                       zeta = angle - (tan(5*PI/12*180/PI) - 5*PI/12),
                       zetavec = [1, zeta, zeta^2, zeta^3, zeta^4, zeta^5, zeta^6],
                       coef = [5*PI/12, 7 - 4*sqrt(3), - (388 - 224*sqrt(3)), 41784 - 24124*sqrt(3),-(5623484 - 3246720*sqrt(3)), 847344640 - 489214656*sqrt(3),
                                - (410335269916/3 - 236907178544*sqrt(3)/3)]
                  )
                  coef*zetavec
            : let(
                   angvec = [1, 1/angle, 1/angle^2, 1/angle^3, 1/angle^4, 1/angle^5],
                   coef = [PI/2, -1, PI/2, -(PI^2/4+2/3), PI^3/8+PI, -(PI^4/16+PI^2+13/15)]
              )
              coef*angvec,
        f = function (inv) let(tan_inv=tan(inv*180/PI)) inv+(angle - (tan_inv - inv))/(tan_inv * tan_inv)
    )
    f(f(f(f(inv))));

Error slightly better than root_find().