"""Harmonic Pancake RevB 30:1 - tooth-skip analysis. Run from anywhere: python analysis.py Reads the three CSVs one folder up, writes figures + tables next to this file. Gear (from gearmaker link): pancake, cycloid profile, m=0.55, z_flex=60, z_circ=62 -> ratio 60/(62-60) = 30. With the output held, one tooth of slip at either mesh costs exactly 1/(z_circ - z_flex) = 0.5 input turns. """ import os import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt HERE = os.path.dirname(os.path.abspath(__file__)) DATA = os.path.dirname(HERE) ARM_M = 0.100 RATIO = 30 TURNS_PER_TOOTH = 0.5 # input turns per tooth slip with output held RUNS = { "baseline": "leftbeef-20261011-122428-207003.csv", # NEMA17 direct drive, no gearbox "run1": "rightbeef-20261011-180808-209421.csv", # 18:08, healthy -> broke "run2": "rightbeef-20261011-205830-100749.csv", # 20:58, same damaged box } # ---- style (dark, reel-friendly) ------------------------------------------- # PLOT_THEME=light renders on the website's cream; PLOT_OUT redirects the figures # (tables and numbers.txt always stay next to this file); PLOT_DPI sets resolution. THEME = os.environ.get("PLOT_THEME", "dark") OUT = os.environ.get("PLOT_OUT", HERE) DPI = int(os.environ.get("PLOT_DPI", "150")) if THEME == "light": SURF, INK, INK2, GRID, TILE = "#faf4e8", "#1c1814", "#4a423a", "#e3d7c2", "#f1e6d2" C_CH1, C_CH2, C_IN, C_SKIP, C_ALT = "#2a78d6", "#eb6834", "#008300", "#e34948", "#4a3aa7" else: SURF, INK, INK2, GRID, TILE = "#1a1a19", "#ffffff", "#c3c2b7", "#3a3a37", "#242422" C_CH1, C_CH2, C_IN, C_SKIP, C_ALT = "#3987e5", "#d95926", "#199e70", "#e66767", "#9085e9" plt.rcParams.update({ "figure.facecolor": SURF, "axes.facecolor": SURF, "savefig.facecolor": SURF, "text.color": INK, "axes.labelcolor": INK2, "xtick.color": INK2, "ytick.color": INK2, "axes.edgecolor": GRID, "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.6, "axes.spines.top": False, "axes.spines.right": False, "lines.linewidth": 2, "font.size": 12, "axes.titlesize": 15, "axes.titleweight": "bold", "axes.titlelocation": "left", "legend.frameon": False, }) def load(key): d = pd.read_csv(os.path.join(DATA, RUNS[key])) d["in_turns"] = np.unwrap(np.radians(d.input_deg)) / (2 * np.pi) d["mt_turns"] = np.unwrap(np.radians(d.output_deg)) / (2 * np.pi) d["T_nm"] = (d.ch1_n - d.ch2_n) * ARM_M # + = pushing CH1, - = pushing CH2 return d def contacts(d, thr=0.3, min_s=0.3): """Contiguous load-cell contacts (per channel) with plateau stats.""" rows = [] t = d.time_s.values for ch in ("ch1_n", "ch2_n"): on = d[ch].values > thr e = np.diff(on.astype(int)) st, en = list(np.where(e == 1)[0] + 1), list(np.where(e == -1)[0] + 1) if on[0]: st.insert(0, 0) if on[-1]: en.append(len(on)) for a, b in zip(st, en): if t[b - 1] - t[a] < min_s: continue mid = slice(a + (b - a) // 3, b - (b - a) // 6) rows.append(dict(ch=ch[:3], t0=t[a], t1=t[b - 1], F_peak=d[ch].values[a:b].max(), F_hold=np.median(d[ch].values[mid]), iq=np.median(np.abs(d.iq_a.values[mid])), cmd=np.median(d.current_cmd_a.values[mid]), in_pos=np.median(d.in_turns.values[mid]))) return pd.DataFrame(rows).sort_values("t0").reset_index(drop=True) def skips(d, f_min=2.5, drop_frac=0.4, advance=0.3, merge_s=0.06): """Tooth-skip events. A skip = the stop force collapses (> drop_frac of its peak within 15 ms) while the motor keeps advancing INTO the stop (> `advance` turns over the next 200 ms). That rules out normal stroke reversals, where the force also drops but the motor backs away. Drops within `merge_s` of a skip are arm bounce, not a new tooth. """ T = d.T_nm.values F = np.abs(T) / ARM_M t, inp = d.time_s.values, d.in_turns.values n = len(F) ev, i = [], 0 while i < n - 15: if F[i] > f_min and F[i + 3] < (1 - drop_frac) * F[i]: j = max(0, i - 10) + np.argmax(F[max(0, i - 10):i + 1]) # true local peak push = np.sign(T[j]) # +: driving into CH1 (input increasing), -: into CH2 if (inp[min(n - 1, j + 40)] - inp[j]) * push > advance: if not ev or t[j] - t[ev[-1]] > merge_s: ev.append(j) i += 4 else: i += 1 # Static skips: arm parked on a stop (motor stalled for 0.2 s), then the motor # lurches > `advance` turns into the stop within 200 ms. The force sags rather # than snaps here, so the drop test above misses it (run 2, 49.5 s). for j in range(40, n - 40): if F[j] < 1.0 or abs(inp[j] - inp[j - 40]) > 0.02: continue push = np.sign(T[j]) if (inp[j + 40] - inp[j]) * push > advance and all(abs(t[j] - t[e]) > 0.3 for e in ev): ev.append(j) ev = np.array(sorted(ev), dtype=int) out = pd.DataFrame(dict(t=t[ev], dir=np.where(T[ev] > 0, "ch1", "ch2"), T_peak=np.abs(T[ev]), in_pos=inp[ev], cmd=d.current_cmd_a.values[ev], iq=d.iq_a.values[ev])) out["d_in"] = out.in_pos.diff().abs() out["main"] = True return out def main(): D = {k: load(k) for k in RUNS} C = {k: contacts(D[k]) for k in RUNS} S = {k: skips(D[k]) for k in ("run1", "run2")} # ---- motor torque constant from direct-drive baseline -------------------- cb = C["baseline"] kt_motor = np.sum(cb.F_hold * ARM_M * cb.iq) / np.sum(cb.iq ** 2) # LSQ through origin ideal = kt_motor * RATIO lines = [] P = lines.append P(f"Motor torque constant (direct drive, LSQ): {kt_motor:.3f} N*m/A -> ideal 30:1 output {ideal:.2f} N*m/A") for k in ("run1", "run2"): c = C[k].copy() c["T_hold"] = c.F_hold * ARM_M c["T_peak"] = c.F_peak * ARM_M c["Nm_per_A"] = c.T_hold / c.iq.replace(0, np.nan) c["static_eff"] = c.Nm_per_A / ideal C[k] = c P(f"\n[{k}] contacts\n" + c.round(3).to_string()) s = S[k] m = s[s.main] P(f"\n[{k}] skip events: {len(s)} total, {len(m)} main (>0.6 N*m)") if len(m) > 1: dd = m.d_in[(m.d_in > 0.3) & (m.d_in < 0.7)] P(f" spacing between consecutive skips: median {dd.median():.3f} input turns " f"(IQR {dd.quantile(.25):.3f}-{dd.quantile(.75):.3f}) -> 1 tooth = {TURNS_PER_TOOTH} turns") P(f" peak torque at skip: median {m.T_peak.median():.2f}, max {m.T_peak.max():.2f}, min {m.T_peak.min():.2f} N*m") P(f" first skip at t={s.t.iloc[0]:.2f}s, T={s.T_peak.iloc[0]:.2f} N*m") C[k].to_csv(os.path.join(HERE, f"contacts_{k}.csv"), index=False) s.to_csv(os.path.join(HERE, f"skips_{k}.csv"), index=False) # position lost = drift of the stop positions in input turns for k in ("run1", "run2"): c = C[k] for ch in ("ch1", "ch2"): p = c[c.ch == ch].in_pos.values P(f"[{k}] {ch} stop input position over time: {np.round(p, 2)}") report = "\n".join(lines) print(report) open(os.path.join(HERE, "numbers.txt"), "w").write(report) d1, s1 = D["run1"], S["run1"] # ---- FIG 1: run1 overview ------------------------------------------------- fig, ax = plt.subplots(2, 1, figsize=(13, 7.5), sharex=True, gridspec_kw=dict(height_ratios=[3, 2])) ax[0].plot(d1.time_s, d1.ch1_n * ARM_M, color=C_CH1, lw=1.2, label="CH1 side") ax[0].plot(d1.time_s, -d1.ch2_n * ARM_M, color=C_CH2, lw=1.2, label="CH2 side") m = s1[s1.main] ax[0].scatter(m.t, np.sign(D["run1"].T_nm.values[np.searchsorted(d1.time_s, m.t)]) * m.T_peak, s=18, color=C_SKIP, zorder=5, label=f"tooth skip ({len(m)})") ax[0].set_ylabel("output torque [N·m]") ax[0].set_title("Run 1 (18:08): the gearbox holds 1.47 N·m, then starts ratcheting") ax[0].legend(loc="upper left", ncol=3) ax[1].plot(d1.time_s, d1.in_turns, color=C_IN, lw=1.5) ax[1].set_ylabel("motor position [turns]") ax[1].set_xlabel("time [s]") for k0, k1 in ((41.4, 43.4), (51.4, 53.4), (56.3, 58.1), (61.4, 63.3)): for a in ax: a.axvspan(k0, k1, color=C_SKIP, alpha=0.12, lw=0) ax[1].annotate(f"stop positions drift:\n{len(s1)} one-tooth skips\n(0.5 motor turns each)", xy=(53.3, -18), xytext=(28, -14), color=INK, fontsize=11, arrowprops=dict(arrowstyle="->", color=INK2)) ax[1].set_xlim(14, 66) fig.tight_layout() fig.savefig(os.path.join(OUT, "fig1_run1_overview.png"), dpi=DPI) # ---- FIG 2: ratchet zoom --------------------------------------------------- z = d1[(d1.time_s > 41.25) & (d1.time_s < 42.35)] fig, ax = plt.subplots(2, 1, figsize=(13, 7.5), sharex=True) ax[0].plot(z.time_s, -z.T_nm, color=C_CH2, lw=2) zs = s1[(s1.t > 41.25) & (s1.t < 42.35) & s1.main] ax[0].scatter(zs.t, zs.T_peak, s=50, color=C_SKIP, zorder=5, label="tooth jumps") ax[0].set_ylabel("torque on CH2 stop [N·m]") ax[0].set_title("Ratcheting: wind up, jump a tooth, repeat, every 0.1 s") ax[0].legend(loc="lower right") ax[1].plot(z.time_s, z.in_turns - z.in_turns.iloc[0], color=C_IN, lw=2, label="motor (input)") ax[1].plot(z.time_s, (z.mt_turns - z.mt_turns.iloc[0]) / 12 * RATIO, color=C_ALT, lw=2, label="output arm (×30, input-equivalent)") for tt, y in zip(zs.t, zs.in_pos - z.in_turns.iloc[0]): ax[1].axvline(tt, color=C_SKIP, lw=0.8, alpha=0.6) ax[1].set_ylabel("position [input turns]") ax[1].set_xlabel("time [s]") ax[1].legend(loc="lower left") ax[1].text(41.27, -1.9, "motor keeps turning…\n…arm doesn't move\n\n0.50 turns between jumps\n= exactly 1 tooth", color=INK, fontsize=12, va="top", bbox=dict(facecolor=SURF, edgecolor="none", pad=4)) fig.tight_layout() fig.savefig(os.path.join(OUT, "fig2_ratchet_zoom.png"), dpi=DPI) # ---- FIG 3: skip spacing histogram ---------------------------------------- allm = pd.concat([S["run1"][S["run1"].main]]) dd = allm.d_in[(allm.d_in > 0.2) & (allm.d_in < 1.2)] fig, ax = plt.subplots(figsize=(10, 5.5)) ax.hist(dd, bins=np.arange(0.2, 1.2, 0.02), color=C_SKIP, edgecolor=SURF, linewidth=1) ax.axvline(0.5, color=INK, ls="--", lw=1.2) ax.text(0.515, ax.get_ylim()[1] * 0.85, "1 tooth = 1/(62−60) = 0.5 motor turns", color=INK) ax.set_xlabel("motor travel between consecutive skips [turns]") ax.set_ylabel("count") ax.set_title(f"Every skip is exactly one tooth (n={len(dd)}, median {dd.median():.3f})") fig.tight_layout() fig.savefig(os.path.join(OUT, "fig3_skip_spacing.png"), dpi=DPI) # ---- FIG 4: capacity ladder ------------------------------------------------ c1 = C["run1"] healthy = c1[(c1.t1 < s1.t.iloc[0]) & (c1.cmd > 0)] s2 = S["run2"] first = s1.iloc[0] later = s1.iloc[1:] ladder = [ ("Run 1: held, no skip", healthy.T_peak.max(), C_IN), (f"Run 1: FIRST skip ({first.t:.1f} s)", first.T_peak, C_SKIP), ("Run 1: ratcheting after that (median)", later.T_peak.median(), C_SKIP), ("Run 2: skips on impact (median)", s2[s2.T_peak > 0.5].T_peak.median(), C_ALT), ("Run 2: skip while parked on stop", s2[s2.T_peak <= 0.5].T_peak.max(), C_ALT), ] fig, ax = plt.subplots(figsize=(12, 5.5)) for i, (lab, v, col) in enumerate(ladder): ax.barh(i, v, color=col, height=0.6) ax.text(v + 0.02, i, f"{v:.2f} N·m", va="center", color=INK, fontsize=13, fontweight="bold") ax.set_yticks(range(len(ladder))) ax.set_yticklabels([l[0] for l in ladder], fontsize=12, color=INK) ax.invert_yaxis() ax.set_xlim(0, 1.75) ax.set_xlabel("output torque [N·m] (= load-cell N × 0.1 m)") ax.grid(axis="y", visible=False) ax.set_title("Once it skips, it keeps getting weaker") fig.tight_layout() fig.savefig(os.path.join(OUT, "fig4_capacity_ladder.png"), dpi=DPI) # ---- FIG 5: torque per amp / static efficiency ----------------------------- fig, ax = plt.subplots(figsize=(10, 6)) ii = np.linspace(0, 1.3, 10) ax.plot(ii, ideal * ii, color=INK2, ls="--", lw=1.5, label=f"ideal 30:1 ({ideal:.1f} N·m/A)") for k, col, lab in (("run1", C_CH1, "run 1 (before skipping)"), ("run2", C_ALT, "run 2 (damaged)")): c = C[k] c = c[(c.cmd > 0)] if k == "run1": c = c[c.t1 < 41.4] ax.scatter(c.iq, c.T_hold, s=70, color=col, label=lab, zorder=4) ax.scatter(cb.iq, cb.F_hold * ARM_M * RATIO, s=40, facecolors="none", edgecolors=C_IN, label="motor alone ×30 (direct-drive test)") ax.set_xlabel("motor current |iq| [A]") ax.set_ylabel("held output torque [N·m]") ax.set_title("Output torque per amp: about half of ideal while stalled") ax.legend(loc="upper left") ax.set_xlim(0, 1.3) ax.set_ylim(0, 5.5) fig.tight_layout() fig.savefig(os.path.join(OUT, "fig5_torque_per_amp.png"), dpi=DPI) # ---- FIG 6: one ratchet cycle, torque vs output-equivalent windup ----------- fig, ax = plt.subplots(figsize=(10, 6)) z = d1[(d1.time_s > 41.40) & (d1.time_s < 42.25)] ax.plot(-(z.in_turns - z.in_turns.iloc[0]) * 360 / RATIO, -z.T_nm, color=C_CH2, lw=1.8) ax.set_xlabel("motor travel, output-equivalent [deg] (6° = 1 tooth)") ax.set_ylabel("torque on stop [N·m]") ax.set_title("Every 6° of motor travel: wind up to ~1.2 N·m, jump one tooth") for k in range(0, 9): ax.axvline(k * 6, color=GRID, lw=1) fig.tight_layout() fig.savefig(os.path.join(OUT, "fig6_sawtooth_vs_angle.png"), dpi=DPI) print("\nfigures written to", HERE) if __name__ == "__main__": main()