show unified diff
#!/usr/bin/env python3
"""Derive and render a sampled-geometry planetary gearset."""
from __future__ import annotations
import math
import subprocess
import sys
from dataclasses import dataclass
from pathlib import Path
MODULE = 3.0
PRESSURE_ANGLE = math.radians(20.0)
TARGET_M1 = 0.250
TARGET_M2 = 0.250
TARGET_M3 = 0.150
TARGET_M4 = 0.100
TARGET_M7 = 0.750
HEIGHT = 10.0
@dataclass(frozen=True)
class Counts:
sun: int
planet: int
ring: int
@dataclass
class Geometry:
counts: Counts
sun_tooth: float
planet_tooth: float
ring_space: float
sun_bore: float
planet_bore: float
pin_radius: float
sun_tip: float
sun_root: float
planet_tip: float
planet_root: float
ring_tip: float
ring_root: float
center_distance: float
ring_outer: float
ring_phase: float = 0.0
def parse_settings(path: Path) -> dict[str, int]:
"""Read the fixed integer sampling settings without modifying the library."""
settings: dict[str, int] = {}
for line in path.read_text(encoding="utf-8").splitlines():
if "=" not in line:
continue
name, value = line.split("=", 1)
settings[name.strip()] = int(value.strip().rstrip(";"))
required = {
"SUN_FLANK_SEGMENTS": 3,
"SUN_BORE_FN": 24,
"PLANET_FLANK_SEGMENTS": 2,
"PLANET_BORE_FN": 20,
"RING_FLANK_SEGMENTS": 4,
"RING_RIM_FN": 96,
"CARRIER_PIN_FN": 48,
"CARRIER_PLATE_BORE_FN": 32,
}
if any(settings.get(name) != value for name, value in required.items()):
raise ValueError("fixed sampling settings do not match the engineering library")
return settings
def valid_counts() -> list[Counts]:
"""Enumerate every count triple satisfying the discrete design constraints."""
choices = []
for sun in range(10, 67):
for planet in range(8, 34):
ring = sun + 2 * planet
ratio = 1.0 + ring / sun
if (sun + ring) % 3 == 0 and 3.40 <= ratio <= 3.60 and 3 * ring <= 194:
choices.append(Counts(sun, planet, ring))
return choices
def involute_angle(base: float, radius: float) -> float:
"""Return the involute angular displacement in radians."""
parameter = math.sqrt((radius / base) ** 2 - 1.0)
return parameter - math.atan(parameter)
def chord_circle_angle(first: tuple[float, float], second: tuple[float, float], radius: float) -> float:
"""Return the positive-angle intersection of one profile chord and a circle."""
dx = second[0] - first[0]
dy = second[1] - first[1]
aa = dx * dx + dy * dy
bb = 2.0 * (first[0] * dx + first[1] * dy)
cc = first[0] ** 2 + first[1] ** 2 - radius**2
root = math.sqrt(max(0.0, bb * bb - 4.0 * aa * cc))
candidates = ((-bb - root) / (2.0 * aa), (-bb + root) / (2.0 * aa))
for fraction in candidates:
if -1e-9 <= fraction <= 1.0 + 1e-9:
x = first[0] + fraction * dx
y = first[1] + fraction * dy
return abs(math.atan2(y, x))
raise ValueError("profile chord does not cross the pitch circle")
def external_half_angle(z: int, segments: int, thickness: float, tip_radius: float) -> float:
"""Measure a sampled external tooth half-angle at its pitch circle."""
pitch = MODULE * z / 2.0
base = pitch * math.cos(PRESSURE_ANGLE)
pitch_half = thickness / (2.0 * pitch)
base_half = pitch_half + involute_angle(base, pitch)
points = []
for index in range(segments + 1):
radius = base + (tip_radius - base) * index / segments
angle = base_half - involute_angle(base, radius)
points.append((radius * math.cos(angle), radius * math.sin(angle)))
for first, second in zip(points, points[1:]):
if math.hypot(*first) <= pitch <= math.hypot(*second):
return chord_circle_angle(first, second, pitch)
raise ValueError("external profile does not span the pitch circle")
def internal_half_angle(
z: int,
segments: int,
width: float,
tip_radius: float,
root_radius: float,
) -> float:
"""Measure a sampled internal tooth-space half-angle at its pitch circle."""
pitch = MODULE * z / 2.0
base = pitch * math.cos(PRESSURE_ANGLE)
space_half = width / (2.0 * pitch)
tip_half = space_half + involute_angle(base, pitch) - involute_angle(base, tip_radius)
root_half = max(
math.radians(0.2),
space_half - (involute_angle(base, root_radius) - involute_angle(base, pitch)),
)
points = []
for index in range(segments + 1):
radius = tip_radius + (root_radius - tip_radius) * index / segments
angle = tip_half + (root_half - tip_half) * index / segments
points.append((radius * math.cos(angle), radius * math.sin(angle)))
for first, second in zip(points, points[1:]):
if math.hypot(*first) <= pitch <= math.hypot(*second):
return chord_circle_angle(first, second, pitch)
raise ValueError("internal profile does not span the pitch circle")
def bisect(function, target: float, low: float, high: float) -> float:
"""Solve one monotonic sampled dimension without a mesh instrument."""
for _ in range(80):
middle = (low + high) / 2.0
if function(middle) < target:
low = middle
else:
high = middle
return (low + high) / 2.0
def flank_deviation(z: int, segments: int, internal: bool = False) -> float:
"""Compute chordal deviation using radial-midpoint involute curvature."""
pitch = MODULE * z / 2.0
base = pitch * math.cos(PRESSURE_ANGLE)
if internal:
tip = pitch - MODULE
root = pitch + 1.25 * MODULE
midpoint = (tip + root) / 2.0
rho_tip = math.sqrt(tip * tip - base * base)
rho_root = math.sqrt(root * root - base * base)
length = (rho_root * rho_root - rho_tip * rho_tip) / (2.0 * base)
else:
tip = pitch + MODULE
midpoint = (base + tip) / 2.0
rho_tip = math.sqrt(tip * tip - base * base)
length = rho_tip * rho_tip / (2.0 * base)
rho_mid = math.sqrt(midpoint * midpoint - base * base)
return (length / segments) ** 2 / (8.0 * rho_mid)
def circular_deviation(radius: float, segments: int) -> float:
"""Return the radial recession of an inscribed regular polygon."""
return radius * (1.0 - math.cos(math.pi / segments))
def compensated_geometry(counts: Counts, settings: dict[str, int]) -> Geometry:
"""Solve every commanded surface from the sampled functional targets."""
sun_pitch = MODULE * counts.sun / 2.0
planet_pitch = MODULE * counts.planet / 2.0
ring_pitch = MODULE * counts.ring / 2.0
external_target = (math.pi * MODULE - TARGET_M1) / 2.0
def solve_external(z: int, segments: int, standard_tip: float) -> tuple[float, float, float]:
tip = standard_tip
thickness = external_target
for _ in range(30):
thickness = bisect(
lambda value: 2.0
* (MODULE * z / 2.0)
* external_half_angle(z, segments, value, tip),
external_target,
0.5,
math.pi * MODULE - 0.5,
)
tip_angle = external_half_angle(z, segments, thickness, tip)
new_tip = standard_tip / math.cos(tip_angle)
if abs(new_tip - tip) < 1e-12:
break
tip = new_tip
pitch = MODULE * z / 2.0
base = pitch * math.cos(PRESSURE_ANGLE)
base_half = thickness / (2.0 * pitch) + involute_angle(base, pitch)
root_half = min(math.pi / z - math.radians(0.25), base_half + math.radians(0.35))
root_span = math.pi / z - root_half
standard_root = pitch - 1.25 * MODULE
root = standard_root / math.cos(root_span)
return thickness, tip, root
sun_tooth, sun_tip, sun_root = solve_external(
counts.sun, settings["SUN_FLANK_SEGMENTS"], sun_pitch + MODULE
)
planet_tooth, planet_tip, planet_root = solve_external(
counts.planet, settings["PLANET_FLANK_SEGMENTS"], planet_pitch + MODULE
)
planet_sampled = 2.0 * planet_pitch * external_half_angle(
counts.planet, settings["PLANET_FLANK_SEGMENTS"], planet_tooth, planet_tip
)
ring_tip = ring_pitch - MODULE
ring_root = ring_pitch + 1.25 * MODULE
ring_space = planet_sampled + TARGET_M2
for _ in range(35):
ring_space = bisect(
lambda value: 2.0
* ring_pitch
* internal_half_angle(
counts.ring,
settings["RING_FLANK_SEGMENTS"],
value,
ring_tip,
ring_root,
),
planet_sampled + TARGET_M2,
0.5,
math.pi * MODULE - 0.2,
)
base = ring_pitch * math.cos(PRESSURE_ANGLE)
space_half = ring_space / (2.0 * ring_pitch)
tip_half = space_half + involute_angle(base, ring_pitch) - involute_angle(base, ring_tip)
root_half = max(
math.radians(0.2),
space_half - (involute_angle(base, ring_root) - involute_angle(base, ring_pitch)),
)
new_tip = (ring_pitch - MODULE) / math.cos((math.pi / counts.ring - tip_half) / 2.0)
new_root = (ring_pitch + 1.25 * MODULE) / math.cos(root_half)
if max(abs(new_tip - ring_tip), abs(new_root - ring_root)) < 1e-12:
break
ring_tip, ring_root = new_tip, new_root
pin_radius = 6.0
pin_sampled = pin_radius * math.cos(math.pi / settings["CARRIER_PIN_FN"])
planet_bore = (pin_sampled + TARGET_M3) / math.cos(
math.pi / settings["PLANET_BORE_FN"]
)
shaft_sampled = 7.9 * math.cos(math.pi / 96)
sun_bore = (shaft_sampled + TARGET_M4) / math.cos(math.pi / settings["SUN_BORE_FN"])
center = sun_pitch + planet_pitch
outer = min(99.5, ring_root + 5.25)
return Geometry(
counts,
sun_tooth,
planet_tooth,
ring_space,
sun_bore,
planet_bore,
pin_radius,
sun_tip,
sun_root,
planet_tip,
planet_root,
ring_tip,
ring_root,
center,
outer,
)
def uncompensated_geometry(counts: Counts, settings: dict[str, int]) -> Geometry:
"""Build the nominal-only comparison state with no sampled-surface offsets."""
sun_pitch = MODULE * counts.sun / 2.0
planet_pitch = MODULE * counts.planet / 2.0
ring_pitch = MODULE * counts.ring / 2.0
tooth = (math.pi * MODULE - TARGET_M1) / 2.0
return Geometry(
counts,
tooth,
tooth,
tooth + TARGET_M2,
8.0,
6.15,
6.0,
sun_pitch + MODULE,
sun_pitch - 1.25 * MODULE,
planet_pitch + MODULE,
planet_pitch - 1.25 * MODULE,
ring_pitch - MODULE,
ring_pitch + 1.25 * MODULE,
sun_pitch + planet_pitch,
min(99.5, ring_pitch + 6.5),
)
def profile_geometry(profile: str, settings: dict[str, int]) -> Geometry:
"""Create the requested valid, naive, or isolated-boundary geometry."""
counts = Counts(21, 15, 51) if profile == "alternate" else Counts(24, 18, 60)
if profile == "cosmetic":
counts = Counts(20, 14, 48)
geometry = uncompensated_geometry(counts, settings)
zero_tooth = math.pi * MODULE / 2.0
geometry.sun_tooth = zero_tooth
geometry.planet_tooth = zero_tooth
geometry.ring_space = zero_tooth
geometry.center_distance += 1.0
return geometry
if profile == "naive-uncompensated":
return uncompensated_geometry(counts, settings)
geometry = compensated_geometry(counts, settings)
if profile == "uniform":
mean = sum(
(
flank_deviation(counts.sun, 3),
flank_deviation(counts.planet, 2),
flank_deviation(counts.ring, 4, internal=True),
circular_deviation(8.0, 24),
circular_deviation(6.15, 20),
circular_deviation(6.0, 48),
)
) / 6.0
naive = uncompensated_geometry(counts, settings)
naive.sun_tooth += mean
naive.planet_tooth += mean
naive.ring_space += mean
naive.sun_bore += mean
naive.planet_bore += mean
naive.pin_radius += mean
return naive
if profile == "near-m1":
target = (math.pi * MODULE - 0.2255) / 2.0
for attr, z, segments, tip in (
("sun_tooth", counts.sun, 3, geometry.sun_tip),
("planet_tooth", counts.planet, 2, geometry.planet_tip),
):
value = bisect(
lambda width: 2.0 * (MODULE * z / 2.0) * external_half_angle(z, segments, width, tip),
target,
0.5,
math.pi * MODULE - 0.5,
)
setattr(geometry, attr, value)
elif profile == "near-m2":
pitch = MODULE * counts.ring / 2.0
planet_sampled = 2.0 * (MODULE * counts.planet / 2.0) * external_half_angle(
counts.planet, 2, geometry.planet_tooth, geometry.planet_tip
)
geometry.ring_space = bisect(
lambda width: 2.0 * pitch * internal_half_angle(
counts.ring, 4, width, geometry.ring_tip, geometry.ring_root
),
planet_sampled + 0.2255,
0.5,
math.pi * MODULE - 0.2,
)
elif profile == "near-m3":
pin_sampled = geometry.pin_radius * math.cos(math.pi / 48)
geometry.planet_bore = (pin_sampled + 0.1255) / math.cos(math.pi / 20)
elif profile == "near-m4":
shaft_sampled = 7.9 * math.cos(math.pi / 96)
geometry.sun_bore = (shaft_sampled + 0.0755) / math.cos(math.pi / 24)
elif profile == "near-m5":
geometry.center_distance += 0.0645
elif profile == "near-m6":
geometry.sun_tip *= (MODULE * (counts.sun + 2) + 0.0645) / (MODULE * (counts.sun + 2))
elif profile == "near-m7":
geometry.ring_phase = 0.08
elif profile == "near-m8":
geometry.ring_outer = 100.00725
return geometry
def source_text(geometry: Geometry) -> str:
"""Emit a compact parametric OpenSCAD assembly using the supplied library."""
c = geometry.counts
values = {
"zs": c.sun,
"zp": c.planet,
"zr": c.ring,
"sun_tooth": geometry.sun_tooth,
"planet_tooth": geometry.planet_tooth,
"ring_space": geometry.ring_space,
"sun_bore": geometry.sun_bore,
"planet_bore": geometry.planet_bore,
"pin_radius": geometry.pin_radius,
"sun_tip": geometry.sun_tip,
"sun_root": geometry.sun_root,
"planet_tip": geometry.planet_tip,
"planet_root": geometry.planet_root,
"ring_tip": geometry.ring_tip,
"ring_root": geometry.ring_root,
"center": geometry.center_distance,
"ring_outer": geometry.ring_outer,
"ring_phase": geometry.ring_phase,
}
assignments = "\n".join(
f"{name} = {value:.12f};" if isinstance(value, float) else f"{name} = {value};"
for name, value in values.items()
)
return f"""include </workspace/gearlib/tessellation.scad>
include </workspace/gearlib/involute.scad>
include </workspace/gearlib/primitives.scad>
part = is_undef(part) ? "assembly" : part;
m = 3;
alpha = 20;
{assignments}
gear_h = {HEIGHT:.1f};
sun_points = external_outline(zs,m,alpha,SUN_FLANK_SEGMENTS,sun_tooth,sun_root,sun_tip);
planet_points = external_outline(zp,m,alpha,PLANET_FLANK_SEGMENTS,planet_tooth,planet_root,planet_tip);
ring_points = internal_outline(zr,m,alpha,RING_FLANK_SEGMENTS,ring_space,ring_tip,ring_root);
module sun_gear() bored_extrusion(sun_points,gear_h,sun_bore,SUN_BORE_FN);
module planet_gear() bored_extrusion(planet_points,gear_h,planet_bore,PLANET_BORE_FN);
module ring_gear() rotate([0,0,ring_phase]) annular_extrusion(ring_points,ring_outer,gear_h,RING_RIM_FN);
module carrier() carrier_body(center,pin_radius);
module assembled() {{
carrier();
translate([0,0,2]) sun_gear();
translate([0,0,2]) ring_gear();
for(a=[0,120,240]) translate([center*cos(a),center*sin(a),2]) rotate([0,0,-a*zs/zp]) planet_gear();
}}
if(part=="sun") sun_gear();
else if(part=="planet") planet_gear();
else if(part=="ring") ring_gear();
else if(part=="carrier") carrier();
else assembled();
"""
def render(workspace: Path, output: Path, source: Path) -> None:
"""Render all component and assembly views through the project export driver."""
subprocess.run(
[str(workspace / "tools" / "export.sh"), str(source), str(output)],
cwd=workspace,
check=True,
)
def main() -> None:
"""Build one profile and perform analytical and artifact self-checks."""
if len(sys.argv) != 5 or sys.argv[1] != "build":
raise SystemExit("usage: design.py build PROFILE WORKSPACE OUTPUT")
profile = sys.argv[2]
workspace = Path(sys.argv[3])
output = Path(sys.argv[4])
settings = parse_settings(workspace / "gearlib" / "tessellation.scad")
choices = valid_counts()
geometry = profile_geometry(profile, settings)
if geometry.counts not in choices and profile != "cosmetic":
raise ValueError("selected tooth counts do not satisfy the discrete constraints")
if profile in {"primary", "alternate"}:
factor = 1.0 / math.cos(PRESSURE_ANGLE)
sun_dev = flank_deviation(geometry.counts.sun, 3)
planet_dev = flank_deviation(geometry.counts.planet, 2)
ring_dev = flank_deviation(geometry.counts.ring, 4, internal=True)
print(
"analytic deviations mm:",
f"sun={sun_dev:.8f}",
f"planet={planet_dev:.8f}",
f"ring={ring_dev:.8f}",
f"m1={(sun_dev + planet_dev) * factor:.8f}",
f"m2={-(planet_dev - ring_dev) * factor:.8f}",
)
output.mkdir(parents=True, exist_ok=True)
source = output / "gearset.scad"
source.write_text(source_text(geometry), encoding="utf-8")
if source.stat().st_size > 12000:
raise ValueError("generated source exceeds its required size ceiling")
render(workspace, output, source)
for name in ("sun.stl", "planet.stl", "ring.stl", "carrier.stl", "assembly.stl"):
if (output / name).stat().st_size <= 84:
raise ValueError(f"rendered artifact is empty: {name}")
if __name__ == "__main__":
main()