# (c) 2014 David Douard # (c) 2023 Jonas Bähr # Based on https://github.com/attoparsec/inkscape-extensions.git # Based on gearUtils-03.js by Dr A.R.Collins # http://www.arc.id.au/gearDrawing.html # # Calculation of Bezier coefficients for # Higuchi et al. approximation to an involute. # ref: YNU Digital Eng Lab Memorandum 05-1 # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License (LGPL) # as published by the Free Software Foundation; either version 2 of # the License, or (at your option) any later version. # for detail see the LICENCE text file. # # FCGear is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Library General Public License for more details. # # You should have received a copy of the GNU Library General Public # License along with FCGear; if not, write to the Free Software # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 from math import cos, sin, tan, pi, acos, asin, atan, sqrt, radians from math import comb as binom def CreateExternalGear(w, m, Z, phi, split=True, addCoeff=1.0, dedCoeff=1.25, filletCoeff=0.375, shiftCoeff=0.0): """ Create an external gear w is wire builder object (in which the gear will be constructed) m is the gear's module (pitch diameter divided by the number of teeth) Z is the number of teeth phi is the gear's pressure angle, in degrees addCoeff is the addendum coefficient (addendum normalized by module) dedCoeff is the dedendum coefficient (dedendum normalized by module) filletCoeff is the root fillet radius, normalized by the module. The default of 0.375 matches the hard-coded value (1.5 * 0.25) of the implementation up to v0.20. The ISO Rack specified 0.38, though. shiftCoeff is the profile shift coefficient (profile shift normalized by module) if split is True, each profile of a teeth will consist in 2 Bezier curves of degree 3, otherwise it will be made of one Bezier curve of degree 4 """ _create_involute_profile( wire_builder=w, module=m, number_of_teeth=Z, pressure_angle=radians(phi), split_involute=split, outer_height_coefficient=addCoeff, inner_height_coefficient=dedCoeff, inner_fillet_coefficient=filletCoeff, profile_shift_coefficient=shiftCoeff) def CreateInternalGear(w, m, Z, phi, split=True, addCoeff=0.6, dedCoeff=1.25, filletCoeff=0.375, shiftCoeff=0.0): """ Create an internal gear w is wire builder object (in which the gear will be constructed) m is the gear's module (pitch diameter divided by the number of teeth) Z is the number of teeth phi is the gear's pressure angle, in degrees addCoeff is the addendum coefficient (addendum normalized by module) The default of 0.6 comes from the "Handbook of Gear Design" by Gitin M. Maitra, with the goal to push the addendum circle beyond the base circle to avoid non-involute flanks on the tips. It in turn assumes, however, that the mating pinion uses a larger value of 1.25. And it's only required for a small number of teeth and/or a relatively large mating gear. Anyways, it's kept here as this was the hard-coded value of the implementation up to v0.20. dedCoeff is the dedendum coefficient (dedendum normalized by module) filletCoeff is the root fillet radius, normalized by the module. The default of 0.375 matches the hard-coded value (1.5 * 0.25) of the implementation up to v0.20. The ISO Rack specified 0.38, though. shiftCoeff is the profile shift coefficient (profile shift normalized by module) if split is True, each profile of a teeth will consist in 2 Bezier curves of degree 3, otherwise it will be made of one Bezier curve of degree 4 """ _create_involute_profile( wire_builder=w, module=m, number_of_teeth=Z, pressure_angle=radians(phi), split_involute=split, rotation=pi/Z, # rotate by half a tooth to align the "inner" tooth with the X-axis outer_height_coefficient=dedCoeff, inner_height_coefficient=addCoeff, outer_fillet_coefficient=filletCoeff, profile_shift_coefficient=shiftCoeff) def _create_involute_profile( wire_builder, module, number_of_teeth, pressure_angle=radians(20.0), split_involute=True, rotation=radians(0), outer_height_coefficient=1.0, inner_height_coefficient=1.0, outer_fillet_coefficient=0.0, inner_fillet_coefficient=0.0, profile_shift_coefficient=0.0): """ Create an involute gear profile in the given wire builder This method can be used to create external as well as internal gear and spline profiles. Thus this method does not use the terms "addednum" and "dedendum" or "tip" and "root", but refers to the elements below the reference circle (i.e. towards the center) as "inner", and those above the reference circle (i.e. away from the center) as "outer". For an external gear, outer_height is the addendum, inner_height is the dedendum, and inner_fillet is the root fillet. For an internal gear, inner_height is the addendum, outer_height is the dedendum, and outer_fillet is the root fillet. The "_coefficient" suffix denotes values normalized by the module. """ profile_shift = profile_shift_coefficient * module outer_height = outer_height_coefficient * module + profile_shift # ext: addendum, int: dedednum inner_height = inner_height_coefficient * module - profile_shift # ext: dedendum, int: addednum # ****** calculate radii # All distances from the center of the profile start with "R". # As those are mostly used in mathematical formulae, we use "symbols" rather than "speaking # names" to keep them short. Otherwise the math becomes unreadable. Rref = number_of_teeth * module / 2 # reference circle radius Rb = Rref * cos(pressure_angle) # base circle radius Ro = Rref + outer_height # radius of outer circle (tip for external, root for internal gears) Ri = Rref - inner_height # radius of inner circle (root for external, tip for internal gears) fi = inner_fillet_coefficient * module # fillet radius at inner circle (ext: root fillet) Rci = Ri + fi # radius at which the center of the inner fillet circle sits Rfi = Rci # radius at which the inner fillet ends (for now assuming some non-involute flank) fo = outer_fillet_coefficient * module # fillet radius at outer circle (int: root fillet) Rco = Ro - fo # radius at which the center of the outer fillet circle sits Rfo = Ro # radius at which the outer fillet ends (for now assuming no outer fillet) has_non_involute_flank = Rfi < Rb # above the base circle we have the involute, # below we have a straight line towards the center. has_inner_fillet = fi > 0 has_outer_fillet = fo > 0 if has_inner_fillet and not has_non_involute_flank: # In this case we need tangency of the involute's start and the inner fillet circle. # It has to be somewhere between max(Rb,Ri) and Rci. # So we need the radius R from the origin depending on the tangent angle, for both: # - the involute: Rinv(t) = Rb * sqrt(1 + t**2), rationale: see below. # - the fillet circle: Rcirc(t) = ... a bit more complex, see below. # and then find the R where both tangents are equal. # As the tangent angle of the involute equals the parameter t, we can just take the # parametric polar representation of the involute as our first equation. # For a circle in parameter form, the tangent angle is pi/2 - t (for the first quadrant, # i.e. 0 <= t <= pi/2, and more than one quadrant is not of interest for us). Unfortunately, # the fillet circle is offset from the origin by (Rci,phi), so things becomes more messy: # Rcirc(t) = sqrt((Rci*cos(phi) + fi*cos(t))**2 + (Rci*sin(phi) + fi*sin(t))**2) # To get rid of the sin(t) part we assume phi "very small", i.e. sin(phi) becomes 0. # This is justyfied because Rci is much larger than fi and the parallax error is # neglectable. Next we substitute the parameter t by pi/2-q-pi. For one to account for the # tangent angle definitions (see above), and then to turn our cicle by 180° as for the # inner fillet we need the third quadrant of the circle (NB: we are looking at the upper # half of the right most tooth, i.e. the involute grows downwards!). This results in # sqrt(2*fi*Rci*cos(-1/2*pi - q) + fi**2 + Rci**2) which simplifies to # Rcirc(q) = sqrt(-2*fi*Rci*sin(q) + fi**2 + Rci**2) # which is the polar r for a given tangent angle q (with the simplification of phi being 0) # This can now be inverted to give the tangent angle at a given r: # Qcirc(r) = asin((-r**2 + fi**2 + Rci**2)/(2*fi*Rci)) # For the involute we have Rinv(q) = Rb * sqrt(1 + q**2), with q=t as the circle is already # adapted, thus: Qinv(r) = +- sqrt(r**2 - Rb**2)/Rb # Now to find the r where the tangents match, our Rfi, we set Qcirc=Qinv and solve for r. # This doesn't seem to have an analytical solution, though, so let's apply Newton's # Method to find the root of q := Qinv - Qcirc over Rb...Rci, or for larger number # of teeth, Ri...Rci as in this case the base circle is below the inner circle. # The graph of q is strictly monotnonous and very steep, no surprises to expect. # To simplify the above equations and to find the derivative, SageMath was used. q = lambda r: (sqrt(r**2 - Rb**2) / Rb - asin((-r**2 + fi**2 + Rci**2) / (2 * fi * Rci))) q_prime = lambda r: (r / (sqrt(-Rb**2 + r**2) * Rb) + r / (fi * Rci * sqrt(1 - 1/4 * (r**2 - fi**2 - Rci**2)**2 / (fi**2 * Rci**2)))) Rfi = findRootNewton(q, q_prime, x_min=max(Rb, Ri), x_max=Rci) if has_outer_fillet: # In this case we need tangency of the involute's end and the outer fillet circle. # For a detailed explanation refer to the inner fillet above. # Here, however, we substitute t with pi/2-q, i.e. we omit the extra -pi, as we are now # interested in the first quadrant: Again, our involute grows downwards but now the fillet # circle needs to approach the involute "from inside the tooth" and at the far end of the # involute. The fillet is now # Rcirc(q) = sqrt(2*fo*Rco*sin(q) + fo**2 + Rco**2) # which can be inversed to # Qcirc(r) = asin((r**2 - fo**2 - Rco**2)/(2*fo*Rco)) # The involute is the very same as for the inner fillet. # However, the simplification of assuming phi=0 is more questionable here as phi doesn't # have an upper bound from the fillet's radius any more but gets larger as the involute # angle grows. Having a non-zero phi in the original Rcirc(q) makes solving impossible, # so lets apply another trick: # If we cannot apply a polar angle to the position of the circle, we could do it for the # involute. Such a rotation doesn't have any influence on Rinv, so not on Qinv, but changes # the interpretation of it: The tangent angle of the involute experiences the same rotation # as the involute itself. So it is just a simple offset: # Our q(r) becomes Qinv(r) - Qcirc(r) - phi_corr, where phi_corr is the amount we # (virtually) rotate our involute around the origin. To bring the fillet circle back on # X-axis, we assume the (real) polar angle of its center position being fo below the # involute's end at Ro. phi_corr = genInvolutePolar(Rb, Ro) + atan(fo / Ro) q = lambda r: (sqrt(r**2 - Rb**2) / Rb - asin((r**2 - fo**2 - Rco**2) / (2 * fo * Rco)) - phi_corr) q_prime = lambda r: (r / (sqrt(-Rb**2 + r**2) * Rb) - r / (fo * Rco * sqrt(1 - 1/4 * (r**2 - fo**2 - Rco**2)**2 / (fo**2 * Rco**2)))) Rfo = findRootNewton(q, q_prime, x_min=max(Rb, Rco), x_max=Ro) # ****** calculate angles (all in radians) angular_pitch = 2 * pi / number_of_teeth # angle subtended by complete tooth/space pair base_to_ref = genInvolutePolar(Rb, Rref) # angle between base and reference circle ref_to_stop = genInvolutePolar(Rb, Rfo) - base_to_ref # angle between ref and involute end if has_non_involute_flank: # involute starts at base circle start_to_ref = base_to_ref else: # involute starts at top of inner fillet, i.e. somewhat above the base circle start_to_ref = base_to_ref - genInvolutePolar(Rb, Rfi) inner_fillet_width = sqrt(fi**2 - (Rci - Rfi)**2) inner_fillet_angle = atan(inner_fillet_width / Rfi) outer_fillet_width = sqrt(fo**2 - (Rfo - Rco)**2) outer_fillet_angle = atan(outer_fillet_width / Rfo) # ****** generate Higuchi involute approximation fe = 1 # fraction of involute length at end of approx fs = 0.01 # fraction of length offset from base to avoid singularity if (not has_non_involute_flank): fs = (Rfi**2 - Rb**2) / (Rfo**2 - Rb**2) # offset start to top of fillet if split_involute: # approximate in 2 sections, split 25% along the involute fm = fs + (fe - fs) / 4 # fraction of length at junction (25% along involute) part1 = BezCoeffs(Rb, Rfo, 3, fs, fm) part2 = BezCoeffs(Rb, Rfo, 3, fm, fe) inv = part1 + part2[1:] # join the 2 sets of coeffs (skip duplicate mid point) else: inv = BezCoeffs(Rb, Rfo, 4, fs, fe) # ****** calculate angular tooth thickness at reference circle enlargement_by_shift = profile_shift * tan(pressure_angle) / Rref tooth_thickness_half_angle = angular_pitch / 4 + enlargement_by_shift psi = tooth_thickness_half_angle # for the formulae below, a symbol is more readable # rotate all points to make the tooth symmetric to the X axis inv = [rotate(pt, -base_to_ref - psi) for pt in inv] # create the back profile of tooth (mirror image on X axis) invR = [mirror(pt) for pt in inv] # ****** calculate section junction points. # Those are the points where the named element ends (start is the end of the previous element). # Suffix _back is back of this tooth, suffix _next is front of next tooth. inner_fillet = toCartesian(Rfi, -psi - start_to_ref) # top of fillet inner_fillet_back = mirror(inner_fillet) # flip to make same point on back of tooth inner_circle_back = toCartesian(Ri, psi + start_to_ref + inner_fillet_angle) inner_circle_next = toCartesian(Ri, angular_pitch - psi - start_to_ref - inner_fillet_angle) inner_fillet_next = rotate(inner_fillet, angular_pitch) # top of fillet, front of next tooth outer_fillet = toCartesian(Ro, -psi + ref_to_stop + outer_fillet_angle) outer_circle = mirror(outer_fillet) # ****** build the gear profile using the provided wire builder thetas = [x * angular_pitch + rotation for x in range(number_of_teeth)] # Make sure we begin *exactly* where our last curve ends. # In theory start == rotate(end, angle_of_last_tooth), but in practice we have limited # precision. Especially if we don't have an inner fillet, we end at inner_circle_next, # not inner_fillet_next. # And even though these two should also be equal (if no inner fillet), they are calculated # differently, which is enough for the resulting wire not being closed any more. # So to be on the save side, we begin at rotate(end, angle_of_last_tooth), not start. if has_inner_fillet: wire_builder.move(rotate(inner_fillet_next, thetas[-1])) else: wire_builder.move(rotate(inner_circle_next, thetas[-1])) for theta in thetas: wire_builder.theta = theta if (has_non_involute_flank): wire_builder.line(inv[0]) # line from inner fillet up to base circle # front involute if split_involute: wire_builder.curve(inv[1], inv[2], inv[3]) wire_builder.curve(inv[4], inv[5], inv[6]) else: wire_builder.curve(*inv[1:]) # is there a section of outer circle between outer fillets? if (outer_circle[1] > outer_fillet[1]): if has_outer_fillet: wire_builder.arc(outer_fillet, fo, 1) # outer fillet wire_builder.arc(outer_circle, Ro, 1) # arc across the outer circle if has_outer_fillet: wire_builder.arc(invR[-1], fo, 1) # outer fillet on back side # back involute if split_involute: wire_builder.curve(invR[5], invR[4], invR[3]) wire_builder.curve(invR[2], invR[1], invR[0]) else: wire_builder.curve(*invR[-2::-1]) if (has_non_involute_flank): wire_builder.line(inner_fillet_back) # line down to top of inner fillet # is there a section of inner circle between inner fillets? if (inner_circle_next[1] > inner_circle_back[1]): if has_inner_fillet: wire_builder.arc(inner_circle_back, fi, 0) # inner fillet on back side wire_builder.arc(inner_circle_next, Ri, 1) # arc across the inner circle if has_inner_fillet: wire_builder.arc(inner_fillet_next, fi, 0) wire_builder.close() def genInvolutePolar(Rb, R): """return the involute angle as function of radius R. Rb = base circle radius """ return (sqrt(R*R - Rb*Rb) / Rb) - acos(Rb / R) def rotate(pt, rads): """rotate pt by rads radians about origin""" sinA = sin(rads) cosA = cos(rads) return (pt[0] * cosA - pt[1] * sinA, pt[0] * sinA + pt[1] * cosA) def mirror(pt): """mirror pt on the X axis, i.e. flip its Y""" return (pt[0], -pt[1]) def toCartesian(radius, angle): """convert polar coords to cartesian""" return (radius * cos(angle), radius * sin(angle)) def findRootNewton(f, f_prime, x_min, x_max): """Apply Newton's Method to find the root of f within x_min and x_max We assume that there is a root in that range and that f is strictly monotonic, i.e. we don't take precautions for overshooting beyond the input range. """ # As initial guess let's take the middle of our input range x = (x_min + x_max) / 2 # FreeCAD.Base.Precision.intersection() is 1e-9, but this file doesn't depend on FreeCAD, # a property very handy when testing in isolation, so let's keep it. PRECISION_INTERSECTION = 1e-9 # Experience has shown that usually 2-3 iterations are enough, but better save than sorry for i in range(6): f_x = f(x) if abs(f_x) < PRECISION_INTERSECTION: return x x = x - f_x/f_prime(x) raise RuntimeError(f"No convergence after {i+1} iterations.") def chebyExpnCoeffs(j, func): N = 50 # a suitably large number N>>p c = 0 for k in range(1, N + 1): c += func(cos(pi * (k - 0.5) / N)) * cos(pi * j * (k - 0.5) / N) return 2 *c / N def chebyPolyCoeffs(p, func): coeffs = [0]*(p+1) fnCoeff = [] T = [coeffs[:] for i in range(p+1)] T[0][0] = 1 T[1][1] = 1 # now generate the Chebyshev polynomial coefficient using # formula T(k+1) = 2xT(k) - T(k-1) which yields # T = [ [ 1, 0, 0, 0, 0, 0], # T0(x) = +1 # [ 0, 1, 0, 0, 0, 0], # T1(x) = 0 +x # [-1, 0, 2, 0, 0, 0], # T2(x) = -1 0 +2xx # [ 0, -3, 0, 4, 0, 0], # T3(x) = 0 -3x 0 +4xxx # [ 1, 0, -8, 0, 8, 0], # T4(x) = +1 0 -8xx 0 +8xxxx # [ 0, 5, 0,-20, 0, 16], # T5(x) = 0 5x 0 -20xxx 0 +16xxxxx # ... ] for k in range(1, p): for j in range(len(T[k]) - 1): T[k + 1][j + 1] = 2 * T[k][j] for j in range(len(T[k - 1])): T[k + 1][j] -= T[k - 1][j] # convert the chebyshev function series into a simple polynomial # and collect like terms, out T polynomial coefficients for k in range(p + 1): fnCoeff.append(chebyExpnCoeffs(k, func)) for k in range(p + 1): for pwr in range(p + 1): coeffs[pwr] += fnCoeff[k] * T[k][pwr] coeffs[0] -= fnCoeff[0] / 2 # fix the 0th coeff return coeffs def bezCoeff(i, p, polyCoeffs): '''generate the polynomial coeffs in one go''' return sum(binom(i, j) * polyCoeffs[j] / binom(p, j) for j in range(i+1)) def BezCoeffs(baseRadius, limitRadius, order, fstart, fstop): """Approximates an involute using a Bezier-curve Parameters: baseRadius - the radius of base circle of the involute. This is where the involute starts, too. limitRadius - the radius of an outer circle, where the involute ends. order - the order of the Bezier curve to be fitted e.g. 3, 4, 5, ... fstart - fraction of distance along the involute to start the approximation. fstop - fraction of distance along the involute to stop the approximation. """ Rb = baseRadius Ra = limitRadius ta = sqrt(Ra * Ra - Rb * Rb) / Rb # involute angle at the limit radius te = sqrt(fstop) * ta # involute angle, theta, at end of approx ts = sqrt(fstart) * ta # involute angle, theta, at start of approx p = order # order of Bezier approximation def involuteXbez(t): "Equation of involute using the Bezier parameter t as variable" # map t (0 <= t <= 1) onto x (where -1 <= x <= 1) x = t * 2 - 1 # map theta (where ts <= theta <= te) from x (-1 <=x <= 1) theta = x * (te - ts) / 2 + (ts + te) / 2 return Rb * (cos(theta) + theta * sin(theta)) def involuteYbez(t): "Equation of involute using the Bezier parameter t as variable" # map t (0 <= t <= 1) onto x (where -1 <= x <= 1) x = t * 2 - 1 # map theta (where ts <= theta <= te) from x (-1 <=x <= 1) theta = x * (te - ts) / 2 + (ts + te) / 2 return Rb * (sin(theta) - theta * cos(theta)) # calc Bezier coeffs bzCoeffs = [] polyCoeffsX = chebyPolyCoeffs(p, involuteXbez) polyCoeffsY = chebyPolyCoeffs(p, involuteYbez) for i in range(p + 1): bx = bezCoeff(i, p, polyCoeffsX) by = bezCoeff(i, p, polyCoeffsY) bzCoeffs.append((bx, by)) return bzCoeffs