Source code for madcad.gear

# This file is part of pymadcad,  distributed under license LGPL v3
'''
	This module provide functions to generate racks and gears with many different shapes.

	The mainstream gears are involute gears (using involute of circles as contact curves). This is due to the standard shape of racks. And involute gears are the best choice for most designs, which is why the current implementation focusses on them.

	However the tools here are modular, which means you can combine them with your own functions to create customized gears.

	Examples
	--------

	As a first use, you may want to generate a fully finished spur gear
	If you do not specify optional arguments, the function will provide good defaults for it.

		>>> # fully featured gear
		>>> gear(step=3, z=12, depth=4, bore_radius=1.5)

	You may want to generate a more bizarre gear, with a a different hub for instance. You will then need to assemble a gear from its 3 components:  exterior (tooths), structure (between hub and exterior), and hub (what holds the gear on an axis)

		>>> # this assemble a gear specifying each independant sub-parts
		>>> ext_radius = (3 * 12) / (2 * pi) - 3
		>>> int_radius = 4
		>>> geargather(
		...	 gearexterior(repeat_circular(gearprofile(3, 30), 30), depth=4),
		...	 gearstructure("rounded", ext_radius, int_radius, 2, patterns=6),
		...	 my_hub_mesh,
		... )

	For reuse in your custom functions, the functions used to generate the gears are exposed:

		>>> # this is the raw profile of a tooth
		>>> gearprofile(step=3, z=12)
'''

from .mathutils import *
from .mesh import Web, Wire, Mesh, web, wire
from .blending import junction, blendpair
from .generation import extrusion, revolution, repeat, extrans
from .primitives import Circle, ArcCentered, Segment
from .triangulation import triangulation
from .selection import *
from .rendering import show
from .cut import bevel, chamfer
from .boolean import intersection
from .io import cachefunc
from . import settings

from math import *
from functools import reduce, partial
from operator import add
from copy import copy
from madcad.mathutils import COMPREC


[docs]def rackprofile(step, h=None, offset=0, alpha=radians(20), resolution=None) -> Wire: ''' Generate a 1-period tooth profile for a rack Parameters: step: period length over the primitive line h: tooth half height offset: rack reference line offset with the primitive line (as a distance) - the primitive line is the adherence line with gears - the reference line is the line half the tooth is above and half below alpha: angle of the tooth sides, a.k.a pressure angle of the contact ''' if h is None: h = default_height(step, alpha) e = offset # change name for convenience x = 0.5 + 2*e/step*tan(alpha) # fraction of the tooth above the primitive circle return Wire([ vec3(step*x/2 - tan(alpha) * ( h-e), h-e, 0), vec3(step*x/2 - tan(alpha) * (-h-e), -h-e, 0), vec3(step*(2-x)/2 + tan(alpha) * (-h-e), -h-e, 0), vec3(step*(2-x)/2 + tan(alpha) * ( h-e), h-e, 0), vec3(step*(2+x)/2 - tan(alpha) * ( h-e), h-e, 0), ], groups=['rack'])
[docs]def gearprofile(step, z, h=None, offset=0, alpha=radians(20), resolution=None, **kwargs) -> Wire: ''' Generate a 1-period tooth profile for a straight gear Parameters: step: period length over the primitive circle z: number of tooth on the gear this profile is meant for h: tooth half height offset: offset distance of the matching rack profile (see above) alpha: pressure angle in the contact ''' if h is None: h = default_height(step, alpha) p = step*z / (2*pi) # primitive circle c = p * cos(alpha) # tangent circle to gliding axis e = offset # change name for convenience #e = offset + 0.05*h # the real offset is 5% (the 1*m, 1.25*m) when the standard says it is 0 x = 0.5 + 2*offset/step*tan(alpha) # fraction of the tooth above the primitive circle o0 = angle(involute(c, 0, tan(alpha))) # offset of contact curve oi = atan((h-e)/p * tan(alpha)) # offset of interference curve l0 = involuteat(c, p+h+e) # interval size of contact curve # if the tooth height is unreachable, find the intersection between the two contact curves t0 = -step/(2*p)*x - o0 t = t0+l0 if involute(c, t0, t)[1] > 0: # Newton solver method for i in range(5): f = sin(t) + cos(t)*(t0-t) J = - sin(t)*(t0-t) t -= f/J l0 = t-t0 # if there is an interference curve interference = c > p-h+e if interference: # Newton solver method # to compute the parameters (t1, t2) of the intersection between contact line and interference line t0, ti = o0, oi # initial state t1 = t0 - tan(alpha) # put contact line on primitive t2 = ti + sqrt(c**2 - (p-h+e)**2) /p # put interference point on base circle for i in range(8): ct1, ct2 = cos(t1), cos(t2) st1, st2 = sin(t1), sin(t2) # function value f = ( c*vec2(ct1,st1) + c*(t0-t1)*vec2(-st1,ct1) + (h-e-p)*vec2(ct2,st2) -p*(ti-t2)*vec2(-st2,ct2) ) # jacobian matrix (f partial derivatives) J = mat2( -c*(t0-t1)*vec2(ct1,st1), p*(ti-t2)*vec2(ct2,st2) + (h-e)*vec2(-st2,ct2), ) # iteration t1, t2 = vec2(t1,t2) - inverse(J)*f li = t2 - ti # interval size of interference curve s0 = t0 - t1 # generation start of contact curve else: s0 = involuteat(c, p-h+e) points = [] tracks = [] n = 2 + settings.curve_resolution(h, step/p, resolution) # number of points to place # parameter for first side place = step/(2*p)*x # place of intersection with the primitive circle t0 = place + o0 # start of contact curve ti = place + oi # start of interference curve # contact line for i in range(n+1): t = interpol1(t0-l0, t0-s0, i/n) v = involute(c, t0, t) points.append(vec3(v,0)) tracks.append(0) # interference line if interference: for i in range(n+1): t = interpol1(ti+li, ti, i/n) v = involuteof(p, ti, -h+e, t) points.append(vec3(v,0)) tracks.append(1) tracks[-1] = 2 # parameters for second side place = step/(2*p)*(2-x) t0 = place - o0 ti = place - oi # interference line if interference: for i in range(n+1): t = interpol1(ti, ti-li, i/n) v = involuteof(p, ti, -h+e, t) points.append(vec3(v,0)) tracks.append(3) # contact line for i in range(n+1): t = interpol1(t0+s0, t0+l0, i/n) v = involute(c, t0, t) points.append(vec3(v,0)) tracks.append(4) tracks[-1] = 5 points.append(angleAxis(step/p, vec3(0,0,1)) * points[0]) tracks.append(5) return Wire(points, tracks=tracks, groups=['contact', 'interference', 'bottom', 'interference', 'contact', 'top'])
[docs]def gearcircles(step, z, h=None, offset=0, alpha=radians(30)): ''' return the convenient circles radius for a gear with the given parameters return is `(primitive, base, bottom, top)` ''' if h is None: h = default_height(step, alpha) e = h*e p = step*z / (2*pi) # primitive circle c = p * cos(alpha) # tangent circle to gliding axis return p, c, p-h-e, p+h-e
def default_height(step, alpha): return 0.5 * 2.25 * step/pi * cos(alpha)/cos(radians(20))
[docs]def involute(c, t0, t): ''' give a point of parameter `t` on involute from circle or radius `c`, starting from `t0` on the circle `t` and `t0` are angular positions ''' x = vec2(cos(t), sin(t)) y = vec2(-x[1], x[0]) return c*(x + y*(t0-t))
[docs]def involuteat(c, r): ''' give the parameter for the involute of circle radius `c` to reach radius `r` ''' return sqrt((r/c)**2 - 1)
[docs]def involuteof(c, t0, d, t): ''' give a point of parameter `t` on involute with offset, from circle or radius `c`, starting from `t0` on the circle `t` and `t0` are angular positions ''' x = vec2(cos(t), sin(t)) y = vec2(-x[1], x[0]) return (c+d)*x + c*y*(t0-t)
def angle(p): return atan2(p[1], p[0])
[docs]def repeat_circular(profile, n: int) -> Wire: ''' Repeat n times the given Wire by rotation around (O,Z) ''' result = repeat(profile, n, rotatearound( anglebt(noproject(profile[0],Z), noproject(profile[-1],Z)), (O,Z), )) result.mergeclose() return result
[docs]def pattern_full(ext_radius, int_radius, depth, int_height=0, **kwargs) -> Mesh: """ Generate two full parts of the structure (the top and the bottom). Return a tuple (Web, Web, None) where the first `Web` is the top of the structure and the second `Web` is the bottom of the structure. Parameters: ext_radius (float): float (radius of the external border of the pattern int_radius (float): radius of the internal border of the pattern depth (float): face width int_height (float): if you want a pinion with a structure thinner than the value of `depth`, the total height will be `total_height = depth - 2 * int_height` """ half_depth = depth / 2 assert half_depth > int_height, "`int_height` must be smaller than `depth / 2`" axis = (O, Z) circles_ref = web(Circle(axis, ext_radius)) + web(Circle(axis, int_radius)).flip() top_profile = circles_ref.transform(vec3(0, 0, half_depth - int_height)) bottom_profile = circles_ref.transform(vec3(0, 0, -half_depth + int_height)).flip() if ext_radius == int_radius: return top_profile + bottom_profile mesh = triangulation(top_profile) + triangulation(bottom_profile) mesh.mergeclose() return mesh
[docs]def pattern_circle( ext_radius, int_radius, depth, int_height = 0, ratio = 1, patterns: int = 5, circles_radius = None, circles_place = None, **kwargs, ) -> Mesh: """ Generate two parts of the structure (the top and the bottom) with `patterns` distributed on the whole structure. Return a tuple (Web, Web, Mesh) where the first `Web` is the top of the structure, the second `Web` is the bottom of the structure and the last element `Mesh` is all side surfaces. Parameters: ext_radius (float): radius of the external border of the structure int_radius (float): radius of the internal border of the structure depth (float): face width int_height (float): if you want a pinion with a structure thinner than the value of `depth`, the total height will be `total_height = depth - 2 * int_height` ratio (float): number that allows to change proportionally the radius of circles patterns (int): number of circles of the structure circles_radius (float): radius of circles circles_place (float): radius where the origins of circles are placed Note: - For instance, with a ratio of 1.5, the radius of circles `circles_radius` will be divided by 1.5 - If `circles_radius` is chosen, `ratio` won't impact the radius of circles `circles_radius` """ # Parameters half_depth = depth / 2 if not (circles_radius) and not (circles_place): assert 0.5 < ratio, "`ratio` must be in the interval ]0.5; +inf[" if not circles_place: circles_place = (ext_radius + int_radius) / 2 if not circles_radius: circles_radius = (ext_radius - int_radius) / (4 * ratio) assert circles_radius / circles_place < sin(pi / patterns), "too many circles for the specified size. Change either 'ratio', 'circle_radius' or 'patterns'" Z = vec3(0,0,1) axis = (vec3(0,0,0), Z) # Borders ext_circle = web(Circle(axis, ext_radius)) int_circle = web(Circle(axis, int_radius)).flip() # Pattern circle_ref = web(Circle((vec3(circles_place, 0, 0), Z), circles_radius)) angle = 2 * pi / patterns pattern_profile = repeat(circle_ref.flip(), patterns, rotatearound(angle, axis)) # Profiles and surfaces webs_ref = (pattern_profile, ext_circle, int_circle) top_webs = [wire.transform(vec3(0, 0, half_depth - int_height)) for wire in webs_ref] bottom_webs = [wire.transform(vec3(0, 0, -half_depth + int_height)).flip() for wire in webs_ref] top_profile = reduce(add, top_webs) bottom_profile = reduce(add, bottom_webs) surfaces = extrusion(vec3(0, 0, depth - 2 * int_height), bottom_webs[0].flip()) mesh = triangulation(top_profile) + triangulation(bottom_profile) + surfaces mesh.mergeclose() return mesh
def create_pattern_rect( ext_radius, int_radius, depth, int_height = 0, ratio = 1, patterns: int = 5, ext_thickness = None, int_thickness = None, rounded: bool = False, **kwargs, ) -> Mesh: """ Function used for `pattern_rect` when `pattern` = "rect" and `pattern_rounded` when `pattern` = "rounded" Parameters: ext_radius (float): radius of the external border of the structure int_radius (float): radius of the internal border of the structure depth (float): face width int_height (float): if you want a pinion with a structure thinner than the value of `depth`, the total height will be `total_height = depth - 2 * int_height` ratio (float): it is a number which is the proportional factor for an homothety of patterns patterns (int): number of patterns inside the structure ext_thickness (float): internal radius of the pattern int_thickness (float): external radius of the pattern rounded (bool): if it is True, the pattern will be rounded Note: - For instance, if `ratio` is 1.5, the area of the pattern will be divided by 1.5. - If `r_int` and `r_ext` are chosen, `ratio` won't impact these parameters """ # Parameters half_depth = depth / 2 offset = ratio * 0.05 * (ext_radius - int_radius) r_ext = ext_radius - (ext_thickness or offset) r_int = int_radius + (int_thickness or offset) thickness = ratio * 0.15 * r_ext # thickness of "arms" # angle between the bottom edge of an arm and the middle of it theta_1 = asin(thickness / (2 * r_int)) # angle between the top edge of an arm and the middle of it theta_2 = asin(thickness / (2 * r_ext)) angle_step = 2 * pi / patterns O = vec3(0) Z = vec3(0, 0, 1) axis = (O, Z) # Borders ext_circle = web(Circle(axis, ext_radius)) int_circle = web(Circle(axis, int_radius)).flip() #function to convert cylindrical coordinates to cartesian coordinates cyl2cart = lambda r, theta: vec3(r * cos(theta), r * sin(theta), 0) # Vertex points of the pattern A1 = cyl2cart(r_int, theta_1) A2 = cyl2cart(r_ext, theta_2) B1 = cyl2cart(r_int, angle_step - theta_1) B2 = cyl2cart(r_ext, angle_step - theta_2) # Primitives pattern_ref = web([ Segment(B2, B1), ArcCentered((O, -Z), B1, A1), Segment(A1, A2), ArcCentered((O, Z), A2, B2), ]) # Multiply the pattern of reference if rounded: # rounded case bevel( pattern_ref, pattern_ref.frontiers(0,1,2), ('radius', min(0.3*(r_ext - r_int), r_int*(angle_step/2-theta_1))), ) bevel( pattern_ref, pattern_ref.frontiers(2,3,0), ('radius', min(0.3*(r_ext - r_int), r_ext*(angle_step/2-theta_2))), ) pattern_ref.mergeclose() pattern_profile = repeat(web(pattern_ref), patterns, rotatearound(angle_step, axis)) # Profiles and surfaces webs_ref = (pattern_profile.flip(), ext_circle, int_circle) top_webs = [web.transform(vec3(0, 0, half_depth - int_height)) for web in webs_ref] bottom_webs = [web.transform(vec3(0, 0, -half_depth + int_height)).flip() for web in webs_ref] top_profile = reduce(add, top_webs) bottom_profile = reduce(add, bottom_webs) surfaces = extrusion(vec3(0, 0, depth - 2 * int_height), bottom_webs[0].flip()) mesh = triangulation(top_profile) + triangulation(bottom_profile) + surfaces mesh.mergeclose() return mesh
[docs]def pattern_rect( ext_radius, int_radius, depth, int_height = 0, ratio = 1, patterns = 5, ext_thickness = None, int_thickness = None, **kwargs, ) -> Mesh: """ Generate two parts of the structure (the top and the bottom) with `patterns` distributed on the whole structure. All corners are straight. Check the function `pattern_rounded` to get rounded corners. Return a tuple (Web, Web, Mesh) where the first `Web` is the top of the structure, the second `Web` is the bottom of the structure and the last element `Mesh` is all side surfaces. Parameters: ext_radius: float (radius of the external border of the structure) int_radius: float (radius of the internal border of the structure) depth: float (face width) int_height: float (if you want a pinion with a structure thinner than the value of `depth`, the total height will be `total_height = depth - 2 * int_height`) ratio: float (it is a number which is the proportional factor for an homothety of patterns) patterns: int (number of patterns inside the structure) ext_thickness (float): internal radius of the pattern int_thickness (float): external radius of the pattern Note: - For instance, if `ratio` is 1.5, the area of the pattern will be divided by 1.5. - If `r_int` and `r_ext` are chosen, `ratio` won't impact these parameters """ return create_pattern_rect(ext_radius, int_radius, depth, int_height, ratio, patterns, ext_thickness, int_thickness, False)
[docs]def pattern_rounded( ext_radius, int_radius, depth, int_height = 0, ratio = 1, patterns: int = 5, ext_thickness = None, int_thickness = None, **kwargs, ) -> Mesh: """ Generate two parts of the structure (the top and the bottom) with `patterns` distributed on the whole structure. All corners are rounded. Check the function `pattern_rect` to get straight corners. Return a tuple (Web, Web, Mesh) where the first `Web` is the top of the structure, the second `Web` is the bottom of the structure and the last element `Mesh` is all side surfaces. Parameters: ext_radius: float (radius of the external border of the structure) int_radius: float (radius of the internal border of the structure) depth: float (face width) int_height: float (if you want a pinion with a structure thinner than the value of `depth`, the total height will be `total_height = depth - 2 * int_height`) ratio: float (it is a number which is the proportional factor for an homothety of patterns) patterns: int (number of patterns inside the structure) ext_thickness (float): internal radius of the pattern int_thickness (float): external radius of the pattern Note: - For instance, if `ratio` is 1.5, the area of the pattern will be divided by 1.5. - If `r_int` and `r_ext` are chosen, `ratio` won't impact these parameters """ return create_pattern_rect(ext_radius, int_radius, depth, int_height, ratio, patterns, ext_thickness, int_thickness, True)
def minmax_radius(points): min, max = inf, 0 for p in points: r = length2(vec2(p)) if r > max: max = r elif r < min: min = r return sqrt(min), sqrt(max)
[docs]def gearexterior( profile: Wire, z, depth, step=None, helix_angle=0, chamfer=0, resolution=None, **kwargs, ) -> Mesh: """ Generate the external part of the pinion Parameters: profile (Web): profile of the pinion generated by `gearprofile` depth (float): extrusion eight - width of the gear along its axis step (float): step of chordal pitch, must be specified for non-null helix_angle, unused otherwise helix_angle (float): helix angle for helical gears - only without bevel; `bevel` = False - it must be a radian angle chamfer (bool | float | (float, float)): set the parameter of chamfer `(angle, ratio)` such as `angle` is the chamfer angle, `ratio` is where the chamfer is applied (`rmin + ratio * (rmax - rmin)`) """ # Parameters half_depth = depth / 2 footer, header = minmax_radius(profile.points) R = header # the largest radius of profile # The code generates a tooth and a body, applies # a boolean operation and then repeat it to # get a complete gear surface # Body refpoint = profile[0] - half_depth * Z # Offsets needed for boolean operation more_offset = (header + 0.3 * (header - footer)) / header less_offset = (footer - 0.3 * (header - footer)) / header A = (more_offset * (X + Y) + Z) * refpoint # external radius offset B = (less_offset * (X + Y) + Z) * refpoint # internal radius offset bottom = web(Segment(A, B)) if helix_angle: # helical tooth case # Step to generate a helical transformation step = settings.curve_resolution( depth / cos(helix_angle), abs(depth * tan(helix_angle) / R), resolution ) angle = depth * tan(helix_angle) / R / (step + 1) h = depth / (step + 1) bottom_gear_edge = profile.transform(translate(-half_depth * Z)) transformations = ( translate(i * h * Z) * rotate(angle * i, Z) for i in range(step + 2) ) links = ((i, i + 1, 0) for i in range(step + 1)) # Generation of helical gear surface bottom_gear_edge = profile.transform(translate(-half_depth * Z)) gear_surface = extrans(bottom_gear_edge, transformations, links) gear_surface.mergeclose() # Transformation from step 0 to step -1 tra = translate(depth * Z) * rotate(depth * tan(helix_angle) / R, Z) # Generation of the body if chamfer: if isinstance(chamfer, (tuple, list)): chamfer_angle, ratio = chamfer elif isinstance(chamfer, float): chamfer_angle, ratio = chamfer, 0.5 else: chamfer_angle, ratio = pi / 8, 0.5 assert chamfer_angle > 0, "chamfer angle must be positive." # Transformation for one step from step 0 to step 1 trs = translate(h * Z) * rotate(angle, Z) # Axis for rotation to generate the chamfer axis = normalize(cross(Z, A)) coeff = (footer + (1 - ratio) * (header - footer)) / header M = (coeff * (X + Y) + Z) * refpoint # Compute closed points C = rotatearound(-chamfer_angle, (M, axis)) * A # top D = rotatearound(chamfer_angle, (M, axis)) * A # bottom # Adjust points t = _get_intersection(refpoint, trs * refpoint, M, D) I = refpoint + t * (trs * refpoint - refpoint) D = M + normalize(I - M) * length(M - D) I = refpoint - t * (trs * refpoint - refpoint) C = M + normalize(I - M) * length(M - C) bottom = web([Segment(C, M), Segment(M, B)]) bottom.mergeclose() top = web([Segment(D, M), Segment(M, B)]).transform(tra) top.mergeclose() else: top = bottom.transform(tra) body = revolution(2 * pi / z, (O, Z), top + bottom.flip()) body.mergeclose() else: # straight tooth case # Generation of body if chamfer: if isinstance(chamfer, (tuple, list)): chamfer_angle, ratio = chamfer elif isinstance(chamfer, float): chamfer_angle, ratio = chamfer, 0.5 else: chamfer_angle, ratio = pi / 8, 0.5 assert chamfer_angle > 0, "chamfer angle must be positive." # Axis for rotation to generate the chamfer axis = normalize(cross(Z, A)) coeff = (footer + (1 - ratio) * (header - footer)) / header M = (coeff * (X + Y) + Z) * refpoint # Compute closed points C = rotatearound(-chamfer_angle, (M, axis)) * A # top D = rotatearound(chamfer_angle, (M, axis)) * A # bottom bottom = web([Segment(C, M), Segment(M, B)]) bottom.mergeclose() top = web([Segment(D, M), Segment(M, B)]).transform(depth * Z) top.mergeclose() else: top = bottom.transform(translate(depth * Z)) body = revolution(2 * pi / z, (O, Z), top + bottom.flip()) body.mergeclose() # Surfaces gear_surface = extrusion( depth * 1.05 * Z, profile.transform(translate(-half_depth * 1.05 * Z)), ) # Boolean operation result = intersection(body.flip(), gear_surface) # Repetition for a complete exterior surface mesh = repeat(result, z, rotatearound(2 * pi / z, (O, Z))) mesh.mergeclose() return mesh
[docs]def gearstructure( pattern, ext_radius, int_radius, depth, int_height = 0, ratio = 1, **kwargs, ) -> Mesh: """ Generate the internal part of the pinion Parameters: ext_radius (float): given by the attribut `_radius` of the result of `repeat_circular` - to avoid interference, it must be smaller than `_radius` (for instance 0.95 * `_radius`)) int_radius (float): it is the same radius of the largest radius of the hub part depth (float): face width pattern: any of 'full', 'circle', 'rect', 'rounded' int_height (float): if you want a pinion with a structure thinner than the value of `depth`, the total height will be `total_height = depth - 2 * int_height` """ # int_radius must not be 0 int_radius = int_radius or 0.1 * ext_radius pattern = "full" if ext_radius == int_radius else pattern pattern_function = globals()['pattern_' + pattern] return pattern_function(ext_radius, int_radius, depth, int_height, ratio=ratio, **kwargs)
[docs]def gearhub( bore_radius, depth, int_height = 0, hub_height = None, hub_radius = None, resolution = None, **kwargs, ) -> Web: """ Generate a hub for a pinion part Parameters: bore_radius (float): radius of the central bore depth (float): face width; same parameter for `gearexterior` and `gearstructure` int_height (float): only useful for no hub case, checkout the function `gearstructure` for more information hub_height (float): height of the hub hub_radius (float): external radius of the hub Note: - if `bore_radius` is null, the function will return a top circle and a bottom circle used for `geargather` function - if `hub_height` is null, the function will return a structure with a bore and without a hub """ if not (bore_radius): # No hub case half_depth = depth / 2 circle_ref = web(Circle((O,Z), depth * 0.1, resolution=resolution)) circle_ref = triangulation(circle_ref) top_surface = circle_ref.transform(vec3(0, 0, half_depth - int_height)) bottom_surface = circle_ref.transform(vec3(0, 0, -half_depth + int_height)).flip() return top_surface + bottom_surface if hub_height is None: hub_height = bore_radius if not (hub_radius): hub_radius = 2 * bore_radius height = hub_height + depth / 2 if hub_height else depth / 2 - int_height axis = (vec3(0, 0, height), vec3(0, 0, 1)) bore_circle = Circle(axis, bore_radius, resolution=resolution) hub_circle = Circle(axis, hub_radius, resolution=resolution) # Top part top_surface = triangulation(web(bore_circle).flip() + web(hub_circle)) # Lateral surfaces lateral_surfaces = extrusion(vec3(0, 0, -hub_height - depth), web(bore_circle)) # Bottom part axis = (vec3(0, 0, -depth / 2), vec3(0, 0, 1)) bottom_bore_circle = Circle(axis, bore_radius, resolution=resolution) bottom_hub_circle = Circle(axis, hub_radius, resolution=resolution) bottom_web = web(bottom_bore_circle) + web(bottom_hub_circle).flip() bottom_surface = triangulation(bottom_web) return top_surface + bottom_surface + lateral_surfaces
[docs]def geargather(exterior, structure, hub) -> Mesh: """ Gather all parts: exterior, structure, and hub You can obtain them via the provided functions, or generate them yourself. """ exterior.mergeclose() structure.mergeclose() hub.mergeclose() # Parameters int_e_radius = minmax_radius(exterior.points)[0] box = structure.box() height_top = box.max.z height_bot = box.min.z int_radius, ext_radius = minmax_radius(structure.points) box = hub.box() height_h_top = box.max.z height_h_bot = box.min.z ext_h_radius = minmax_radius(hub.points)[1] assert ext_h_radius * COMPREC <= int_radius assert ext_radius * COMPREC <= int_e_radius # borderline overlapping check prec2 = (4*NUMPREC*ext_radius)**2 def radius_height(circle): p = circle.points[circle.edges[0][0]] return vec2(length(p - circle.barycenter()), p.z) def samecircle(c1, c2): return distance2(radius_height(c1), radius_height(c2)) < prec2 # select and join all borderlines # j1 is the junction between `exterior` and `structure` # j2 is the junction between `structure` and `hub` top = j1_top = j2_top = bottom = j1_bot = j2_bot = Mesh() circle_int_e_top = select(exterior, vec3(0, 0, height_top)) circle_ext_h_top = select(hub, vec3(ext_h_radius, 0, height_h_top)) if not samecircle(circle_int_e_top, circle_ext_h_top): circle_ext_s_top = select(structure, vec3(ext_radius, 0, height_top)) circle_int_s_top = select(structure, vec3(int_radius, 0, height_top)) if not samecircle(circle_ext_s_top, circle_int_e_top): j1_top = junction(circle_ext_s_top, circle_int_e_top, tangents="straight", resolution=('div',0)) if not samecircle(circle_ext_h_top, circle_int_s_top): j2_top = junction(circle_ext_h_top, circle_int_s_top, tangents="straight", resolution=('div',0)) top = j1_top + j2_top circle_int_e_bot = select(exterior, vec3(0, 0, height_bot)) circle_ext_h_bot = select(hub, vec3(ext_h_radius, 0, height_h_bot)) if not samecircle(circle_int_e_bot, circle_ext_h_bot): circle_ext_s_bot = select(structure, vec3(ext_radius, 0, height_bot)) circle_int_s_bot = select(structure, vec3(int_radius, 0, height_bot)) if not samecircle(circle_ext_s_bot, circle_int_e_bot): j1_bot = junction(circle_ext_s_bot, circle_int_e_bot, tangents="straight", resolution=('div',0)) if not samecircle(circle_ext_h_bot, circle_int_s_bot): j2_bot = junction(circle_ext_h_bot, circle_int_s_bot, tangents="straight", resolution=('div',0)) bottom = j1_bot + j2_bot if isinstance(structure, Web): mesh = exterior + top + bottom + hub else: mesh = exterior + top + structure + bottom + hub return mesh.finish()
[docs]@cachefunc def gear( step, z: int, depth, bore_radius = 0, int_height = 0, pattern = 'full', **kwargs, ) -> Mesh: """ Generate a pinion. Any extra argument will go to functions `gearprofile`, `gearstructure`, or `gearhub` Parameters: step (float): tooth step over the primitive curve, same as the matching rack step the primitive perimeter is `step * z` and radius is `step * z / 2*pi` z (int): number of teeth depth (float): extrusion eight - width of the gear along its axis bore_radius (float): radius of the main bore int_height (float): if you want a pinion with a structure thinner than the value of `depth`, the total height will be `total_height = depth - 2 * int_height` pattern: determine the structure between exterior (tooth) and hub This argument specifies the use a a function named `'pattern_'+pattern` in this module. * Extra parameters for `gearprofile` offset (float): offset of tooth (as a distance) alpha (float): pressure angle in radian * Extra parameters for `gearexterior` helix_angle (float): helix angle to get a helical pinion in radian chamfer (bool | float | (float, float)): set the parameter of chamfer `(angle, ratio)` such as `angle` is the chamfer angle, `ratio` is where the chamfer is applied (`rmin + ratio * (rmax - rmin)`) * Extra parameters for `gearstructure` ratio (float): influence the proportion of dimensions of the structure Note: `int_height` impacts the thickness of the structure unless specified * Extra parameters for `gearhub` hub_height (float): height of the hub shoulder Note: - `int_height` impacts the height of the hub unless specified. - if `hub_height` is null, there will be no hub. """ # profile = repeat_circular(gearprofile(step, z, **kwargs), z) profile = gearprofile(step, z, **kwargs) # Parts exterior = gearexterior(profile, z, depth, step, **kwargs) ext_int = minmax_radius(exterior.points)[0] hub = gearhub(bore_radius, depth, int_height, hub_radius=min(2 * bore_radius, 0.9 * ext_int), **kwargs) hub_ext = minmax_radius(hub.points)[1] structure = gearstructure( pattern, 0.95 * ext_int, max(hub_ext, 0.2 * ext_int), depth, int_height, **kwargs) return geargather(exterior, structure, hub)
def get_pitch_cone_angle(z_pinion:int, z_wheel:int, shaft_angle:float=0.5 * pi) -> float: """ Return the pitch cone angle of the pinion called `gamma_p`. The pitch cone angle of the wheel is equal to `shaft_angle - gamma_p` Parameters: z_pinion (int): the number of teeth on the bevel pinion z_wheel (int): the number of teeth on the bevel wheel shaft_angle (float): the shaft angle """ return atan2(sin(shaft_angle), ((z_wheel / z_pinion) + cos(shaft_angle)))
[docs]def spherical_involute(cone_angle:float, t0:float, t:float) -> vec3: """ Return spherical involute function Parameters: t (float): the angular position t0 (float): the difference phase cone_angle (float): the cone angle Return: a normalized `vec3` """ cos_g, sin_g = cos(cone_angle), sin(cone_angle) return vec3( sin_g * cos(t * sin_g) * cos(t + t0) + sin(t * sin_g) * sin(t + t0), sin_g * cos(t * sin_g) * sin(t + t0) - sin(t * sin_g) * cos(t + t0), cos_g * cos(t * sin_g), )
[docs]def spherical_involuteof(pitch_cone_angle:float, t0:float, alpha:float, t:float) -> vec3: """ Return the spherical interference function Parameters: t (float): the angular position t0 (float): the difference phase pitch_cone_angle (float): the pitch cone angle alpha(float): the height angle offset of the rack Return: a normalized `vec3` """ cos_p, sin_p = cos(pitch_cone_angle), sin(pitch_cone_angle) return ( + cos(alpha) * spherical_involute(pitch_cone_angle, t0, t) + sin(alpha) * vec3(-cos_p * cos(t+t0), -cos_p * sin(t+t0), sin_p) )
def derived_spherical_involute(cone_angle:float, t0:float): """ Return the function of the derived spherical involute function. Parameters: cone_angle (float): the cone angle t0 (float): the phase difference """ cos_g, sin_g = cos(cone_angle), sin(cone_angle) return lambda t: vec3( cos_g ** 2 * sin(t * sin_g) * cos(t + t0), cos_g ** 2 * sin(t * sin_g) * sin(t + t0), -cos_g * sin_g * sin(t * sin_g), ) def jacobian_spherical_involute(base_cona_angle:float, pitch_cone_angle:float, t01:float, t02:float, alpha:float) -> callable: """ Return the function of the jacobian used for the newton method in `spherical_gearprofile` Parameters: base_cona_angle (float): the base cone angle pitch_cone_angle (float): the pitch cone angle t01 (float): the phase of the spherical involute function t02 (float): the phase of the spherical interference function alpha (float): the height angle offset of the rack """ dsi = derived_spherical_involute # for convenience derived_involute = dsi(base_cona_angle, t01) cos_p = cos(pitch_cone_angle) vec = lambda t: vec3(cos_p * sin(t), -cos_p * cos(t), 0) derived_interference = lambda t: dsi(pitch_cone_angle, t02)(t) * cos(alpha) + sin(alpha) * vec(t + t02) return lambda t1, t2: mat3(derived_involute(t1), -derived_interference(t2), vec3(0, 0, 1)) def spherical_rack_tools(z:float, pressure_angle:float=pi / 9, ka:float=1, kd:float=1.25): """ Return a list of all information useful to generate a spherical rack. Five elements : 1) the minimum abscissa for the function (fifth element) 2) the maximum abscissa for the function (fifth element) 3) the phase of a tooth 4) the phase of space 5) the function to generate a tooth Parameters : z (float): number of tooth of the rack equal to `z_pinion / sin(pitch_cone_angle)` or `z_wheel / sin(shaft_angle - pitch_cone_angle)` pressure_angle (float): the pressure angle of the gear ka (float): the addendum coefficient kd (float): the dedendum coefficient """ k = 1 / z gamma_p = 0.5 * pi gamma_b = asin(cos(pressure_angle)) cos_b, sin_b = cos(gamma_b), sin(gamma_b) gamma_f = gamma_p + 2 * ka * k gamma_r = gamma_p - 2 * kd * k phi_p = acos(tan(gamma_b) / tan(gamma_p)) theta_p = atan2(sin_b * tan(phi_p), 1) / sin_b - phi_p phase_diff = k * pi + 2 * theta_p t_min = acos(cos(gamma_r) / cos_b) / sin_b t_max = acos(cos(gamma_f) / cos_b) / sin_b involute = lambda t, t0: spherical_involute(gamma_b, t0, t) v = vec3(1, 1, 0) phase_empty = 2 * pi * k - anglebt(involute(t_min, 0) * v, involute(-t_min, phase_diff) * v) return [t_min, t_max, phase_diff, phase_empty, involute]
[docs]def spherical_rackprofile(z:float, pressure_angle:float=pi / 9, ka:float=1, kd:float=1.25, resolution=None): """ Return a `Wire` which is a tooth of the rack. Parameters: z (float): number of tooth of the rack equal to `z_pinion / sin(pitch_cone_angle)` or `z_wheel / sin(shaft_angle - pitch_cone_angle)` pressure_angle (float): the pressure angle of the gear ka (float): the addendum coefficient kd (float): the dedendum coefficient """ div = settings.curve_resolution(1, pi/z, resolution) t_min, t_max, phase1, phase2, involute = spherical_rack_tools(z, pressure_angle, ka, kd) side1 = [involute(t, 0) for t in linrange(t_min, t_max, div=div)] side2 = [involute(-t, phase1) for t in linrange(t_min, t_max, div=div)] segment = Segment(involute(t_min, -phase2), side1[0]).mesh() segment2 = Segment(side1[-1], side2[-1]).mesh() wire = segment + Wire(side1) + segment2 + Wire(side2).flip() return wire
[docs]def spherical_gearprofile( z:int, pitch_cone_angle:float, pressure_angle:float=pi / 9, ka:float=1, kd:float=1.25, resolution = None, ) -> Wire: """ Generate and return a `Wire` of a 1-period tooth spherical profile for a bevel gear Parameters: z (int): number of tooth on the gear this profile is meant for pitch_cone_angle (float): the pitch cone angle pressure_angle (float): pressure angle of the tooth ka (float): addendum coefficient kd (float): dedendum coefficient """ # Initialization of parameters gamma_p = pitch_cone_angle # for convenience gamma_b = asin(cos(pressure_angle) * sin(gamma_p)) cos_b, sin_b = cos(gamma_b), sin(gamma_b) tooth_size = pi / z involute = lambda t, t0 : spherical_involute(gamma_b, t0, t) epsilon_p = acos(cos(gamma_p) / cos_b) / sin_b theta_p = anglebt(involute(0, 0) * vec3(1, 1, 0), involute(epsilon_p, 0) * vec3(1, 1, 0)) phase_diff = tooth_size + 2 * theta_p phase_empty = phase_interference = 2 * pi / z - phase_diff # The following number `k` is useful to simplify some calculations # It's broadly speaking `1/z_rack` and `z_rack` is not an integer ! k = sin(gamma_p) / z # Spherical involute part gamma_f = gamma_p + 2 * ka * k # addendum cone angle gamma_r = gamma_p - 2 * kd * k # dedendum cone angle t_min = 0 t_max = acos(cos(gamma_f) / cos_b) / sin_b if gamma_r > gamma_b: v = vec3(1, 1, 0) t_min = acos(cos(gamma_r) / cos_b) / sin_b phase_empty = 2 * pi / z - anglebt( involute(t_min, 0) * v, involute(-t_min, phase_diff) * v ) # Calculation of offsets due to geometry of spherical rack _, t_rack_max, phase1, _, rinvolute = spherical_rack_tools(1 / k, pressure_angle, ka, kd) interference = lambda t, t0 : spherical_involuteof(gamma_p, t0, alpha, t) alpha = 2 * ka * k n1, n2 = rinvolute(t_rack_max, 0) * vec3(1, 1, 0), rinvolute(-t_rack_max, phase1) * vec3(1, 1, 0) beta = 0.5 * anglebt(n1, n2) * length(n1) / sin_b # Newton method to calculate the intersection between # the spherical involute and the spherical interference. # Objective function involuteat = lambda t2, t0 : spherical_involuteof(gamma_p, t0, alpha, t2) f = lambda t1, t2: involute(t1, 0) - involuteat(t2, -0.5 * phase_interference + beta) # Jacobian matrix J = jacobian_spherical_involute(gamma_b, gamma_p, 0, -0.5 * phase_interference + beta, alpha) # Compute the intersection values t1, t2, t3 = 0.5 * t_max, -0.5 * t_max, 0 for i in range(8): t1, t2, t3 = vec3(t1, t2, t3) - inverse(J(t1, t2)) * f(t1, t2) # Build sides of a tooth div = settings.curve_resolution(1, 2*pi/z, resolution) interference1 = [interference(t, -0.5 * phase_interference + beta) for t in linrange(0, t2, div=div)] interference2 = [interference(-t, phase_diff + 0.5 * phase_interference - beta) for t in linrange(0, t2, div=div)] side1 = Wire(interference1[:-1]) + Wire([involute(t, 0) for t in linrange(t1, t_max, div=div)]) side2 = Wire(interference2[:-1]) + Wire([involute(-t, phase_diff) for t in linrange(t1, t_max, div=div)]) side2.groups = side1.groups = ['interference', 'contact'] # Extreme points of sides to compute angle between them a = interference(0, -0.5 * phase_interference + beta) b = interference(0, phase_diff + 0.5 * phase_interference - beta) final_phase_empty = 2 * pi / z - anglebt(a * vec3(1, 1, 0), b * vec3(1, 1, 0)) top = Segment(involute(t_max, 0), involute(-t_max, phase_diff)).mesh() bottom = Segment(angleAxis(-final_phase_empty, vec3(0, 0, 1)) * a, a).mesh() top.groups = ['top'] bottom.groups = ['bottom'] return bottom + side1 + top + side2.flip()
def cone_projection(profile: Wire, pitch_cone_angle:float) -> Wire: """ Return a Wire of the spherical profile projected on a cone Parameters: profile (Wire) : the spherical profile pitch_cone_angle (float) : the pitch cone angle """ ref = lambda t: vec3(sin(pitch_cone_angle) * cos(t), sin(pitch_cone_angle) * sin(t), cos(pitch_cone_angle)) new_points = [1 / dot(ref(atan2(point.y, point.x)), point) * point for point in profile.points] return Wire(new_points, indices=profile.indices) # @cachefunc
[docs]def bevelgear( step:float, z:int, pitch_cone_angle:float, pressure_angle:float=pi/9, ka:float=1, kd:float=1.25, helix_angle:float=None, bore_radius:float=None, bore_height:float=None, resolution = None, ): """ Generate a bevel gear. Parameters: step (float): tooth step over the primitive curve, same as the matching rack step the primitive perimeter is `step * z` and radius is `step * z / (2 * pi)` z (int): number of teeth pitch_cone_angle (float): pitch cone angle pressure_angle (float): the pressure angle of the tooth ka (float) : addendum coefficient kd (float) : dedendum coefficient helix_angle (float): helix angle of the tooth bore_radius (float): radius of the main bore bore_height (float): height of the main bore """ if helix_angle: return helical_bevel_gear(step, z, pitch_cone_angle, pressure_angle, ka, kd, helix_angle, bore_radius, bore_height, resolution) else: return straight_bevel_gear(step, z, pitch_cone_angle, pressure_angle, ka, kd, bore_radius, bore_height, resolution)
def straight_bevel_gear( step:float, z:int, pitch_cone_angle:float, pressure_angle:float=pi/9, ka:float=1, kd:float=1.25, bore_radius:float=None, bore_height:float=None, resolution = None, ): """ Generate a bevel gear where teeth are straight. Parameters: step (float): tooth step over the primitive curve, same as the matching rack step the primitive perimeter is `step * z` and radius is `step * z / (2 * pi)` z (int): number of teeth pitch_cone_angle (float): pitch cone angle pressure_angle (flaot): the pressure angle of the tooth ka (float) : addendum coefficient kd (float) : dedendum coefficient bore_radius (float): radius of the main bore bore_height (float): height of the main bore """ # Initialization of parameters gamma_p = pitch_cone_angle # for convenience gamma_b = asin(cos(pressure_angle) * sin(gamma_p)) cos_b, sin_b = cos(gamma_b), sin(gamma_b) rp = z * step / (2 * pi) rho1 = rp / sin(gamma_p) rho0 = 2 * rho1 / 3 k = sin(gamma_p) / z gamma_r = gamma_p - 2 * kd * k gamma_f = gamma_p + 2 * ka * k phi_p = acos(tan(gamma_b) / tan(gamma_p)) theta_p = atan2(sin_b * tan(phi_p), 1) / sin_b - phi_p phase_diff = pi / z + 2 * theta_p # Generate spherical profiles spherical_profile = spherical_gearprofile(z, gamma_p, pressure_angle, ka, kd, resolution=resolution) # one tooth # Generate teeth border teeth_border = extrusion((rho1*1.1)/(rho0*0.9), spherical_profile.transform(rho0*0.9)) # Common values v = vec3(1, 1, 0) # angle1tooth = anglebt(spherical_profile[0] * v, spherical_profile[-1] * v) angle1tooth = 2 * pi / z gamma_l = gamma_p + 2 * (ka + 0.25) * k # offset : 0.25 for boolean operation A = vec3(rho1 * sin(gamma_r), 0, rho1 * cos(gamma_r)) B = vec3(rho1 * sin(gamma_l), 0, rho1 * cos(gamma_l)) C = vec3(rho0 * sin(gamma_r), 0, rho0 * cos(gamma_r)) D = vec3(rho0 * sin(gamma_l), 0, rho0 * cos(gamma_l)) v1 = A * v v2 = spherical_profile[0] * v phase_tooth_body = quat(v1, v2) # Generate points for a section if bore_radius is None: bore_radius = 0.5 * rho0 * sin(gamma_r) if bore_height is None: bore_height = 0.4 * rho1 * cos(gamma_r) if bore_radius: if bore_height: E = vec3(1.5 * bore_radius, 0, rho1 * cos(gamma_r)) F = vec3(bore_radius, 0, rho0 * cos(gamma_r)) G = vec3(1.5 * bore_radius, 0, rho1 * cos(gamma_r) + bore_height) # top of the bore H = vec3(bore_radius, 0, rho1 * cos(gamma_r) + bore_height) # top of the bore wire = Wire([D, C, F, H, G, E, A, B]).segmented() chamfer(wire, [2, 3, 4], ("distance", bore_radius * 0.05)) bevel(wire, [5], ("distance", bore_height * 0.1)) else: E = vec3(bore_radius, 0, rho1 * cos(gamma_r)) F = vec3(bore_radius, 0, rho0 * cos(gamma_r)) wire = Wire([D, C, F, E, A, B]).segmented() chamfer(wire, [2, 3], ("distance", bore_radius * 0.05)) else: E = vec3(0, 0, rho1 * cos(gamma_r)) F = vec3(0, 0, rho0 * cos(gamma_r)) wire = Wire([D, C, F, E, A, B]).segmented() axis = (O, Z) body = revolution(angle1tooth, axis, wire, resolution=resolution) one_tooth = intersection(body.transform(phase_tooth_body), teeth_border) all_teeth = repeat(one_tooth, z, rotatearound(angle1tooth, axis)) all_teeth.finish() # For better orientation : the first tooth is placed centered in the plane (O, Y) # try `show([mybevelgear, square((O, Y), size_of_mybevelgear)])` involute = lambda t, t0: spherical_involute(gamma_b, t0, t) t_max = acos(cos(gamma_f) / cos_b) / sin_b middle_tooth = 0.5 * (involute(t_max, 0) + involute(-t_max, phase_diff)) phase = anglebt(X, middle_tooth * v) return all_teeth.transform(angleAxis(-phase, Z)) def helical_bevel_gear( step:float, z:int, pitch_cone_angle:float, pressure_angle:float = radians(20), ka:float=1, kd:float=1.25, helix_angle:float = radians(20), bore_radius:float=None, bore_height:float=None, resolution = None, ): """ Generate a bevel gear where teeth are helical. Parameters: step (float): tooth step over the primitive curve, same as the matching rack step the primitive perimeter is `step * z` and radius is `step * z / (2 * pi)` z (int): number of teeth pitch_cone_angle (float): pitch cone angle pressure_angle (float): the pressure angle of the tooth ka (float) : addendum coefficient kd (float) : dedendum coefficient helix_angle (float): helix angle of the tooth bore_radius (float): radius of the main bore bore_height (float): height of the main bore """ # Initialization of parameters gamma_p = pitch_cone_angle # for convenience gamma_b = asin(cos(pressure_angle) * sin(gamma_p)) cos_b, sin_b = cos(gamma_b), sin(gamma_b) rp = z * step / (2 * pi) rho1 = rp / sin(gamma_p) rho0 = 2 * rho1 / 3 k = sin(gamma_p) / z gamma_r = gamma_p - 2 * kd * k # dedendum angle gamma_f = gamma_p + 2 * ka * k # addendum angle gamma_l = gamma_p + 2 * (ka + 0.25) * k # offset : 0.25 for boolean operation phi_p = acos(tan(gamma_b) / tan(gamma_p)) theta_p = atan2(sin_b * tan(phi_p), 1) / sin_b - phi_p phase_diff = pi / z + 2 * theta_p angle1tooth = 2 * pi / z v = vec3(1, 1, 0) # useful for computation # Generate spherical profiles tooth = spherical_gearprofile(z, gamma_p, pressure_angle, ka, kd, resolution=resolution) # one tooth # One helical tooth # General parameters cosp = cos(gamma_p) z0 = rho0 * cosp z1 = rho1 * cosp delta = sin(gamma_p) / tan(helix_angle) # parameter of the conical helix angle_step = lambda z: (log(z) - log(z0)) / delta # Compute the number of steps for discretization topbot_angle = (log(z1) - log(z0)) / delta step = settings.curve_resolution( (z1 - z0) / (cos(helix_angle) * cosp), # Length of the conical helix abs(topbot_angle), # Angle of the conical helix resolution, ) # Parameters for scaling and rotation scale_step = (rho1 - rho0) / (step + 1) z_step = (z1 - z0) / (step + 1) # The transformation aims to generate a conical helix. helix_scale = lambda x : scale(vec3(x)) helix_rotate = lambda x: rotate(angle_step(x), Z) gear_surface = extrans(tooth, transformations = ( helix_scale(i * scale_step + rho0) * helix_rotate(z_step * i + z0) for i in range(-1, step + 2 + 1)), links = ( (i, i + 1, 0) for i in range(step + 1 + 2)), ) # Get angle between the initial body and tooth body_point = vec3(sin(gamma_r), 0, cos(gamma_r)) tooth_point = tooth[0] v1 = body_point * v v2 = tooth_point * v phase_tooth_body = anglebt(v1, v2) phase_tooth_body = -phase_tooth_body if cross(v1, v2).z > 0 else phase_tooth_body # Initialization of first points of `body` # A is the first point of the spherical profile with the transformation at step `-1` # B is the first point of spherical profile with the transformation at step `0` # C is the first point of body on the dedendum radius of the spherical profile # D is the first point of body on the adendum radius of the spherical profile # Note: due to discretization, there is an issue of imperfection # A new point E must be computed to solve this problem quat_t2b = angleAxis(-phase_tooth_body, Z) def get_body_points(rho, z, angle): """ Return two points of body according to the radius of sphere and the height. They are useful to build the body correctly for the boolean operation """ mat4_scale_rot = mat4(quat_t2b) * helix_scale(rho) * helix_rotate(z) gear_surface_point_A = helix_scale(rho - scale_step) * helix_rotate(z - z_step) * tooth_point gear_surface_point_B = helix_scale(rho) * transform(angleAxis(angle, Z)) * tooth_point body_C = mat4_scale_rot * vec3(sin(gamma_r), 0, cos(gamma_r)) body_D = mat4_scale_rot * vec3(sin(gamma_l), 0, cos(gamma_l)) # Compute intersection and get E t_rho = _get_intersection(gear_surface_point_A, gear_surface_point_B, body_C, body_D) I = gear_surface_point_A + t_rho * (gear_surface_point_B - gear_surface_point_A) E = body_C + normalize(I - body_C) * length(body_D - body_C) return body_C, E C, D = get_body_points(rho0, z0, 0) A, B = get_body_points(rho1, z1, topbot_angle) # Generate points for a section if bore_radius is None: bore_radius = 0.5 * rho0 * sin(gamma_r) if bore_height is None: bore_height = 0.4 * rho1 * cos(gamma_r) if bore_radius: if bore_height: E = vec3(1.5 * bore_radius, 0, rho1 * cos(gamma_r)) F = vec3(bore_radius, 0, rho0 * cos(gamma_r)) G = vec3(1.5 * bore_radius, 0, rho1 * cos(gamma_r) + bore_height) # top of the bore H = vec3(bore_radius, 0, rho1 * cos(gamma_r) + bore_height) # top of the bore wire = Wire([D, C, F, H, G, E, A, B]).segmented() chamfer(wire, [2, 3, 4], ("distance", bore_radius * 0.05)) bevel(wire, [5], ("distance", bore_height * 0.1)) else: E = vec3(bore_radius, 0, rho1 * cos(gamma_r)) F = vec3(bore_radius, 0, rho0 * cos(gamma_r)) wire = Wire([D, C, F, E, A, B]).segmented() chamfer(wire, [2, 3], ("distance", bore_radius * 0.05)) else: E = vec3(0, 0, rho1 * cos(gamma_r)) F = vec3(0, 0, rho0 * cos(gamma_r)) wire = Wire([D, C, F, E, A, B]).segmented() body = revolution(angle1tooth, (O, Z), wire) # show([gear_surface, body]) onetooth = intersection(gear_surface, body) all_teeth = repeat(onetooth, z, rotatearound(angle1tooth, (O, Z))).finish() # align the mesh on the middle of a tooth involute = lambda t, t0: spherical_involute(gamma_b, t0, t) t_max = acos(cos(gamma_f) / cos_b) / sin_b middle_tooth = 0.5 * (involute(t_max, 0) + involute(-t_max, phase_diff)) return all_teeth.transform(quat(X, middle_tooth * v)) def _get_intersection(A: vec3, B: vec3, C: vec3, D: vec3) -> float: """ Return the intersection parameter between the straight line AB and the straight line CE given: - `xE * xE + yE * yE = xD * xD + yD * yD` - `zE = zD` Parameters: A (vec3): Point B (vec3): Point C (vec3): Point D (vec3): Point Example: ``` t = _get_intersection(A, B, C, D) AB = B - A I = A + t * AB # intersection point ``` """ AB = B - A CD = D - C AC = C - A cA = AB.z / CD.z cB = -AC.z / CD.z m1 = AB.x + cA * C.x m2 = AB.y + cA * C.y p1 = -AC.x + cB * C.x p2 = -AC.y + cB * C.y R2 = D.x * D.x + D.y * D.y a = m1 * m1 + m2 * m2 - cA * cA * R2 b = 2 * (m1 * p1 + m2 * p2 - R2 * cA * cB) c = p1 * p1 + p2 * p2 - cB * cB * R2 delta = b * b - 4 * a * c return (-b + sqrt(delta)) / (2 * a) def matrix4placement(z:int, shaft_angle:float) -> mat4x4: """ Return a matrix of transformation given initial state (Z is the initial axis of revolution). This matrix is helpful when you want to place bevel gears. Parameters: z (int): number of tooth on the gear shaft_angle (float): the shaft angle Example: ```python pinion = bevelgear(step, z_pinion, ...) wheel = bevelgear(step, z_wheel, ...) myparts = Solid(pinion=pinion, ...) matrix = matrix4placement(z_pinion, shaft_angle) show([wheel, myparts.transform(matrix)]) ``` """ return mat4(angleAxis(shaft_angle, Y)) * mat4(angleAxis(pi + pi / z, Z))