Cam and Follower Mechanism Experiments

# Cam Experiment 1: SVAJ and the Fundamental Law of Cam Design
import numpy as np
import matplotlib.pyplot as plt
import os

os.makedirs('plots', exist_ok=True)

# Normalized profiles over a unit rise, x in [0,1]. Return S, V, A, J = dS/dx and derivatives.
def svaj(profile, x):
    tp = 2 * np.pi
    if profile == 'const-velocity':
        return x, np.ones_like(x), np.zeros_like(x), np.zeros_like(x)
    if profile == 'const-accel':
        S = np.where(x <= 0.5, 2 * x**2, 1 - 2 * (1 - x)**2)
        V = np.where(x <= 0.5, 4 * x, 4 * (1 - x))
        A = np.where(x <= 0.5, 4.0, -4.0)
        return S, V, A, np.zeros_like(x)
    if profile == 'simple-harmonic':
        return (0.5 * (1 - np.cos(np.pi * x)), (np.pi / 2) * np.sin(np.pi * x),
                (np.pi**2 / 2) * np.cos(np.pi * x), -(np.pi**3 / 2) * np.sin(np.pi * x))
    if profile == 'cycloidal':
        return (x - np.sin(tp * x) / tp, 1 - np.cos(tp * x),
                tp * np.sin(tp * x), tp**2 * np.cos(tp * x))
    raise ValueError(profile)

# Physical scaling: a rise of lift h over angular width beta (deg) at cam speed rpm.
h, beta_deg, rpm = 20.0, 90.0, 180.0
beta = np.deg2rad(beta_deg)
omega = rpm * 2 * np.pi / 60.0
x = np.linspace(0, 1, 600)

profiles = ['const-velocity', 'const-accel', 'simple-harmonic', 'cycloidal']
colors = {'const-velocity': '#ef4444', 'const-accel': '#f97316',
          'simple-harmonic': '#E9C46A', 'cycloidal': '#2A9D8F'}

fig, axes = plt.subplots(2, 2, figsize=(13, 9))
labels = ['Displacement s (mm)', 'Velocity v (mm/s)', 'Acceleration a (mm/s2)', 'Jerk j (mm/s3)']
for profile in profiles:
    S, V, A, J = svaj(profile, x)
    s = h * S
    v = (h / beta) * V * omega
    a = (h / beta**2) * A * omega**2
    j = (h / beta**3) * J * omega**3
    for ax, y in zip(axes.ravel(), [s, v, a, j]):
        ax.plot(x * beta_deg, y, color=colors[profile], lw=2, label=profile)

for ax, lab in zip(axes.ravel(), labels):
    ax.set_xlabel('Angle through rise (deg)'); ax.set_ylabel(lab)
    ax.grid(True, alpha=0.3); ax.legend(fontsize=8)
axes[0, 0].set_title('Displacement looks fine for all four')
axes[1, 0].set_title('Acceleration: the law-breakers jump at the ends')
plt.tight_layout()
plt.savefig('plots/exp1_svaj.png', dpi=150, bbox_inches='tight')
plt.show()

print("Join condition at the dwell boundary (x=0): the rise must start with v=0 AND a=0,")
print("because the dwell it meets has both zero. A step in either one breaks the law.")
for profile in profiles:
    _, V0, A0, _ = svaj(profile, np.array([0.0]))
    v0 = (h / beta) * V0[0] * omega
    a0 = (h / beta**2) * A0[0] * omega**2
    obeys = abs(v0) < 1e-6 and abs(a0) < 1e-6
    print(f"  {profile:16s}: v(0) = {v0:7.1f} mm/s, a(0) = {a0:9.0f} mm/s2  -> "
          f"{'joins smoothly (OK)' if obeys else 'STEP at the dwell join'}")
print("\nOnly cycloidal starts and ends with both v and a at zero, so it joins the dwells")
print("without a step. Constant velocity steps in v (an infinite acceleration spike);")
print("constant acceleration and simple harmonic step in a (an infinite jerk spike).")

# Cam Experiment 2: The Motion-Profile Shoot-Out
import numpy as np
import matplotlib.pyplot as plt
import os

os.makedirs('plots', exist_ok=True)

# Build normalized S, V(=dS/dx), A(=d2S/dx2), J(=d3S/dx3) on a fine grid for each profile.
N = 4000
x = np.linspace(0, 1, N + 1)
tp = 2 * np.pi

def closed(profile):
    if profile == 'cycloidal':
        return (x - np.sin(tp * x) / tp, 1 - np.cos(tp * x), tp * np.sin(tp * x), tp**2 * np.cos(tp * x))
    if profile == '3-4-5':
        return (10*x**3 - 15*x**4 + 6*x**5, 30*x**2 - 60*x**3 + 30*x**4,
                60*x - 180*x**2 + 120*x**3, 60 - 360*x + 360*x**2)
    if profile == '4-5-6-7':
        return (35*x**4 - 84*x**5 + 70*x**6 - 20*x**7,
                140*x**3 - 420*x**4 + 420*x**5 - 140*x**6,
                420*x**2 - 1680*x**3 + 2100*x**4 - 840*x**5,
                840*x - 5040*x**2 + 8400*x**3 - 4200*x**4)
    raise ValueError(profile)

# Modified trapezoid / modified sine: define the unit-amplitude acceleration shape, then
# integrate twice and normalize so S(1)=1 (exactly how the simulator builds them).
def integral_profile(kind):
    if kind == 'mod-trap':
        A0 = np.where(x < 0.125, np.sin(4*np.pi*x),
             np.where(x < 0.375, 1.0,
             np.where(x < 0.625, np.sin(4*np.pi*(x - 0.25)),
             np.where(x < 0.875, -1.0, np.sin(4*np.pi*(x - 1))))))
    else:  # mod-sine
        A0 = np.where(x < 0.125, np.sin(4*np.pi*x),
             np.where(x < 0.875, np.sin((4*np.pi/3)*(x - 0.125) + np.pi/2),
             np.sin(4*np.pi*(x - 1))))
    V0 = np.concatenate([[0], np.cumsum(0.5*(A0[1:] + A0[:-1]) * np.diff(x))])
    S0 = np.concatenate([[0], np.cumsum(0.5*(V0[1:] + V0[:-1]) * np.diff(x))])
    scale = 1.0 / S0[-1]
    return S0*scale, V0*scale, A0*scale  # A amplitude == Ca

profiles = {}
for p in ['cycloidal', '3-4-5', '4-5-6-7']:
    S, V, A, J = closed(p)
    profiles[p] = dict(Cv=np.max(np.abs(V)), Ca=np.max(np.abs(A)), Cj=np.max(np.abs(J)))
for kind, name in [('mod-trap', 'modified trapezoid'), ('mod-sine', 'modified sine')]:
    S, V, A = integral_profile(kind)
    # jerk = dA/dx via gradient (finite Cj for these profiles)
    J = np.gradient(A, x)
    profiles[name] = dict(Cv=np.max(np.abs(V)), Ca=np.max(np.abs(A)), Cj=np.max(np.abs(J)))

print(f"{'Profile':20s} {'Cv':>7s} {'Ca':>8s} {'Cj':>8s}")
for name, c in profiles.items():
    print(f"{name:20s} {c['Cv']:7.4f} {c['Ca']:8.4f} {c['Cj']:8.2f}")

names = list(profiles.keys())
Ca = [profiles[n]['Ca'] for n in names]
Cj = [profiles[n]['Cj'] for n in names]
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
axes[0].bar(names, Ca, color='#2A9D8F'); axes[0].set_ylabel('Peak accel coefficient Ca')
axes[0].set_title('Lower Ca = smaller inertia force')
axes[1].bar(names, Cj, color='#E9C46A'); axes[1].set_ylabel('Peak jerk coefficient Cj')
axes[1].set_title('Lower Cj = smoother, less vibration')
for ax in axes:
    ax.tick_params(axis='x', rotation=30); ax.grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.savefig('plots/exp2_profiles.png', dpi=150, bbox_inches='tight')
plt.show()

best_a = min(profiles, key=lambda n: profiles[n]['Ca'])
best_j = min(profiles, key=lambda n: profiles[n]['Cj'])
print(f"\nLowest peak acceleration: {best_a}")
print(f"Lowest peak jerk:         {best_j}")

# Cam Experiment 3: Pressure Angle and Cam Sizing
import numpy as np
import matplotlib.pyplot as plt
import os

os.makedirs('plots', exist_ok=True)

# Double-dwell program: dwell, cycloidal rise, dwell, cycloidal fall, dwell (degrees).
h = 25.0           # lift (mm)
rr = 10.0          # roller radius (mm)
prog = [('dwell', 45), ('rise', 90), ('dwell', 90), ('fall', 90), ('dwell', 45)]

def cycloidal(x):
    tp = 2 * np.pi
    return x - np.sin(tp*x)/tp, 1 - np.cos(tp*x)   # S, V=dS/dx

def s_and_v(theta_deg):
    """Return follower displacement s (mm) and v = ds/dtheta (mm/rad) at a cam angle."""
    start, lift = 0.0, 0.0
    n = len(prog)
    for i, (kind, dur) in enumerate(prog):
        last = (i == n - 1)
        if (start <= theta_deg < start + dur) or (last and theta_deg >= start):
            if kind == 'dwell':
                return lift, 0.0
            x = (theta_deg - start) / dur
            S, V = cycloidal(x)
            beta = np.deg2rad(dur)
            sign = 1.0 if kind == 'rise' else -1.0
            return lift + sign * h * S, sign * h * V / beta
        if kind == 'rise':
            lift += h
        elif kind == 'fall':
            lift -= h
        start += dur
    return lift, 0.0

def max_pressure_angle(base_radius, eps=0.0):
    Rp = base_radius + rr
    thetas = np.arange(0, 360, 0.5)
    amax = 0.0
    for th in thetas:
        s, v = s_and_v(th)
        denom = s + np.sqrt(max(Rp**2 - eps**2, 1e-9))
        alpha = np.degrees(np.arctan2(v - eps, denom))
        amax = max(amax, abs(alpha))
    return amax

radii = np.array([20, 25, 30, 35, 40, 45, 50, 55, 60], dtype=float)
amax = np.array([max_pressure_angle(rb) for rb in radii])

plt.figure(figsize=(8, 5))
plt.plot(radii, amax, 'o-', color='#2A9D8F', lw=2, label='Max pressure angle')
plt.axhline(30, color='#ef4444', ls='--', label='30 deg limit')
plt.xlabel('Base-circle radius rb (mm)'); plt.ylabel('Max pressure angle (deg)')
plt.title('A bigger base circle lowers the pressure angle')
plt.legend(); plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('plots/exp3_pressure.png', dpi=150, bbox_inches='tight')
plt.show()

below = radii[amax <= 30]
print("Max pressure angle vs base radius:")
for rb, a in zip(radii, amax):
    print(f"  rb={rb:4.0f} mm -> alpha_max = {a:5.1f} deg  {'OK' if a <= 30 else 'too high'}")
if below.size:
    print(f"Smallest base radius keeping alpha <= 30 deg: about {below[0]:.0f} mm")

# Offset effect at a tight base radius: pressure angle through the RISE only.
rb_tight = 22.0
print(f"\nOffset effect at rb = {rb_tight:.0f} mm (rise-side pressure angle near mid-rise):")
for eps in [0, 5, 10]:
    s, v = s_and_v(90.0)   # mid-rise
    denom = s + np.sqrt(max((rb_tight + rr)**2 - eps**2, 1e-9))
    a = np.degrees(np.arctan2(v - eps, denom))
    print(f"  offset={eps:2d} mm -> alpha = {a:5.1f} deg")
print("Positive offset subtracts from the rise velocity term, easing the rise (it costs you on the fall).")

# Cam Experiment 4: Radius of Curvature and Undercutting
import numpy as np
import matplotlib.pyplot as plt
import os

os.makedirs('plots', exist_ok=True)

# Cycloidal double-dwell. A shorter rise makes the curvature tighter (smaller rho_min).
GOLDEN = 2.5

def cyc(x):
    tp = 2*np.pi
    return x - np.sin(tp*x)/tp, 1 - np.cos(tp*x), tp*np.sin(tp*x)  # S, V, A (d2S/dx2)

def svaj_at(prog, theta_deg, h):
    start, lift = 0.0, 0.0
    n = len(prog)
    for i, (kind, dur) in enumerate(prog):
        last = (i == n - 1)
        if (start <= theta_deg < start + dur) or (last and theta_deg >= start):
            if kind == 'dwell' or dur == 0:
                return lift, 0.0, 0.0
            x = (theta_deg - start) / dur
            S, V, A = cyc(x)
            beta = np.deg2rad(dur)
            sign = 1.0 if kind == 'rise' else -1.0
            return lift + sign*h*S, sign*h*V/beta, sign*h*A/beta**2
        if kind == 'rise': lift += h
        elif kind == 'fall': lift -= h
        start += dur
    return lift, 0.0, 0.0

def rho_min(base_radius, rr, h, rise):
    pad = (360 - 100 - 2*rise) / 2
    prog = [('dwell', pad), ('rise', rise), ('dwell', 100), ('fall', rise), ('dwell', pad)]
    Rp = base_radius + rr
    rmin = np.inf
    for th in np.arange(0, 360, 0.5):
        s, v, a = svaj_at(prog, th, h)
        R = Rp + s
        den = R**2 + 2*v**2 - R*a
        if abs(den) < 1e-6:
            continue
        rho = (R**2 + v**2)**1.5 / den
        if 0 < rho < rmin:
            rmin = rho
    return rmin

# (base radius, roller radius, lift, rise angle) spanning generous to undercut
cases = [(45, 8, 20, 110), (35, 10, 25, 90), (25, 12, 30, 70), (16, 16, 35, 50), (12, 22, 40, 45)]
print(f"{'rb':>4s} {'rr':>4s} {'h':>4s} {'rise':>5s} {'rho_min':>9s} {'ratio':>7s}  verdict")
ratios = []
for rb, rr, h, rise in cases:
    rmin = rho_min(rb, rr, h, rise)
    ratio = rmin / rr
    undercut = rmin < rr
    verdict = 'UNDERCUT' if undercut else ('marginal' if ratio < GOLDEN else 'OK')
    ratios.append(ratio)
    print(f"{rb:4.0f} {rr:4.0f} {h:4.0f} {rise:5.0f} {rmin:9.1f} {ratio:7.2f}  {verdict}")

labels = [f'rb{rb}/rr{rr}' for rb, rr, _, _ in cases]
plt.figure(figsize=(9, 5))
bars = plt.bar(labels, ratios, color=['#2A9D8F' if r >= GOLDEN else ('#E9C46A' if r >= 1 else '#ef4444') for r in ratios])
plt.axhline(1.0, color='#ef4444', ls='--', label='Undercut threshold (ratio = 1)')
plt.axhline(GOLDEN, color='#f59e0b', ls=':', label=f'Golden rule (ratio = {GOLDEN})')
plt.ylabel('rho_min / roller radius'); plt.title('Tight base + big roller -> undercut')
plt.legend(); plt.grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.savefig('plots/exp4_curvature.png', dpi=150, bbox_inches='tight')
plt.show()

print("\nUndercut when rho_min < roller radius (ratio < 1). The golden rule keeps the ratio")
print("at or above about 2.5 so the cam is comfortably manufacturable, not just barely possible.")

# Cam Experiment 5: Follower Dynamics and Jump
import numpy as np
import matplotlib.pyplot as plt
import os

os.makedirs('plots', exist_ok=True)

# Cycloidal double-dwell. Follower: single-DOF spring + mass.
prog = [('dwell', 60), ('rise', 100), ('dwell', 40), ('fall', 100), ('dwell', 60)]
h = 20.0           # mm
m = 0.3            # kg effective mass
k = 25000.0        # N/m spring rate
F0_default = 150.0 # N preload

def cyc(x):
    tp = 2*np.pi
    return x - np.sin(tp*x)/tp, tp*np.sin(tp*x)  # S, A (=d2S/dx2)

def s_and_a(theta_deg):
    start, lift = 0.0, 0.0
    n = len(prog)
    for i, (kind, dur) in enumerate(prog):
        last = (i == n-1)
        if (start <= theta_deg < start+dur) or (last and theta_deg >= start):
            if kind == 'dwell':
                return lift, 0.0
            x = (theta_deg - start)/dur
            S, A = cyc(x)
            beta = np.deg2rad(dur)
            sign = 1.0 if kind == 'rise' else -1.0
            return lift + sign*h*S, sign*h*A/beta**2   # s (mm), a (mm/rad2)
        if kind == 'rise': lift += h
        elif kind == 'fall': lift -= h
        start += dur
    return lift, 0.0

thetas = np.arange(0, 360, 0.5)
s_arr = np.array([s_and_a(t)[0] for t in thetas])
a_arr = np.array([s_and_a(t)[1] for t in thetas])   # mm/rad2

def min_contact_force(rpm, F0):
    omega = rpm*2*np.pi/60.0
    a_time = a_arr * omega**2 * 1e-3        # m/s2
    inertia = m * a_time                    # N
    spring = F0 + k * (s_arr*1e-3)          # N
    return np.min(spring + inertia)

def jump_speed(F0):
    # smallest rpm at which min contact force reaches 0
    for rpm in range(50, 6000, 10):
        if min_contact_force(rpm, F0) <= 0:
            return rpm
    return None

# Contact force vs angle at a few speeds (default preload)
plt.figure(figsize=(9, 5))
for rpm in [800, 1500, 2200]:
    omega = rpm*2*np.pi/60.0
    F = (F0_default + k*(s_arr*1e-3)) + m*(a_arr*omega**2*1e-3)
    plt.plot(thetas, F, lw=2, label=f'{rpm} rpm')
plt.axhline(0, color='#ef4444', ls='--', label='Separation (F = 0)')
plt.xlabel('Cam angle (deg)'); plt.ylabel('Contact force (N)')
plt.title('Raising speed pushes the contact force toward zero')
plt.legend(); plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('plots/exp5_jump.png', dpi=150, bbox_inches='tight')
plt.show()

print("Jump speed vs spring preload:")
for F0 in [40, 80, 150, 250]:
    js = jump_speed(F0)
    print(f"  preload={F0:3.0f} N -> jump speed ~ {js if js else '>6000'} rpm")
print("More preload (or rate) raises the jump speed; a higher-acceleration profile lowers it.")

# Cam Experiment 6: Design a Cam to a Specification
import numpy as np
import os

os.makedirs('plots', exist_ok=True)

# Specification
h = 18.0           # mm lift
rise_deg = 80.0
rpm = 600.0
rr = 8.0           # mm roller
omega = rpm*2*np.pi/60.0
ALPHA_MAX, GOLDEN, m, k, F0 = 30.0, 2.5, 0.4, 20000.0, 120.0

prog = [('dwell', 60), ('rise', rise_deg), ('dwell', 100), ('fall', rise_deg),
        ('dwell', 360 - 60 - 2*rise_deg - 100)]

# Normalized profiles (S, V=dS/dx, A=d2S/dx2)
def profile(name, x):
    tp = 2*np.pi
    if name == 'cycloidal':
        return x - np.sin(tp*x)/tp, 1 - np.cos(tp*x), tp*np.sin(tp*x)
    if name == 'mod-trap':
        # build by integrating the unit acceleration shape, normalized so S(1)=1
        xx = np.linspace(0, 1, 4001)
        A0 = np.where(xx < 0.125, np.sin(4*np.pi*xx),
             np.where(xx < 0.375, 1.0,
             np.where(xx < 0.625, np.sin(4*np.pi*(xx-0.25)),
             np.where(xx < 0.875, -1.0, np.sin(4*np.pi*(xx-1))))))
        V0 = np.concatenate([[0], np.cumsum(0.5*(A0[1:]+A0[:-1])*np.diff(xx))])
        S0 = np.concatenate([[0], np.cumsum(0.5*(V0[1:]+V0[:-1])*np.diff(xx))])
        sc = 1.0/S0[-1]
        return (np.interp(x, xx, S0*sc), np.interp(x, xx, V0*sc), np.interp(x, xx, A0*sc))
    raise ValueError(name)

def svaj_at(name, theta_deg):
    start, lift = 0.0, 0.0
    n = len(prog)
    for i, (kind, dur) in enumerate(prog):
        last = (i == n-1)
        if (start <= theta_deg < start+dur) or (last and theta_deg >= start):
            if kind == 'dwell' or dur == 0:
                return lift, 0.0, 0.0
            x = np.array([(theta_deg - start)/dur])
            S, V, A = profile(name, x)
            beta = np.deg2rad(dur); sign = 1.0 if kind == 'rise' else -1.0
            return lift + sign*h*S[0], sign*h*V[0]/beta, sign*h*A[0]/beta**2
        if kind == 'rise': lift += h
        elif kind == 'fall': lift -= h
        start += dur
    return lift, 0.0, 0.0

def evaluate(name, base_radius):
    Rp = base_radius + rr
    thetas = np.arange(0, 360, 0.5)
    amax, rmin, fmin = 0.0, np.inf, np.inf
    for th in thetas:
        s, v, a = svaj_at(name, th)
        denom = s + Rp
        amax = max(amax, abs(np.degrees(np.arctan2(v, denom))))
        R = Rp + s; den = R**2 + 2*v**2 - R*a
        if abs(den) > 1e-6:
            rho = (R**2 + v**2)**1.5/den
            if 0 < rho < rmin: rmin = rho
        F = (F0 + k*s*1e-3) + m*(a*omega**2*1e-3)
        fmin = min(fmin, F)
    return amax, rmin/rr, fmin

print(f"Spec: h={h} mm, rise={rise_deg} deg, {rpm:.0f} rpm, roller {rr} mm")
print(f"Constraints: alpha<= {ALPHA_MAX} deg, rho ratio>= {GOLDEN}, contact force > 0\n")
for name in ['mod-trap', 'cycloidal']:
    print(f"{name}:")
    for rb in [25, 30, 35, 40]:
        amax, ratio, fmin = evaluate(name, rb)
        ok = (amax <= ALPHA_MAX) and (ratio >= GOLDEN) and (fmin > 0)
        print(f"  rb={rb:2d} mm: alpha_max={amax:4.1f} deg, rho ratio={ratio:4.1f}, "
              f"min contact={fmin:6.0f} N  -> {'MEETS SPEC' if ok else 'no'}")
print("\nPick the smallest base radius that meets every constraint with margin.")