Four-Bar Linkage Experiments

import numpy as np
import matplotlib.pyplot as plt
import csv
import os

# ============================================================
# STEP 1: Load your CSV data from the simulator
# Place your downloaded CSV in the data/ folder, or update the path
# ============================================================
csv_file = 'data/exp1_baseline.csv'
sim_data = None

if os.path.exists(csv_file):
    with open(csv_file) as f:
        reader = csv.DictReader(f)
        sim_data = list(reader)
    print(f"Loaded {len(sim_data)} rows from {csv_file}")
else:
    print(f"No CSV found at {csv_file}. Using analytical data only.")
    print("To use simulator data: download CSV from the simulator or")
    print("record values from the charts into a CSV with these columns:")
    print("Input Angle (deg), Coupler Angle (deg), Follower Angle (deg), ...")

# ============================================================
# STEP 2: Analytical computation (for comparison)
# ============================================================
def solve_fourbar(a, b, c, d, e, theta2_deg):
    t2 = np.radians(theta2_deg)
    Ax, Ay = a*np.cos(t2), a*np.sin(t2)
    dx, dy = d - Ax, e - Ay
    dist = np.sqrt(dx**2 + dy**2)
    if dist > b+c or dist < abs(b-c) or dist < 1e-10:
        return None
    ca = (b**2 - c**2 + dist**2) / (2*dist)
    hsq = b**2 - ca**2
    if hsq < 0: return None
    h = np.sqrt(hsq)
    mx, my = Ax + ca*dx/dist, Ay + ca*dy/dist
    Bx, By = mx + h*dy/dist, my - h*dx/dist
    t3 = np.degrees(np.arctan2(By-Ay, Bx-Ax))
    t4 = np.degrees(np.arctan2(By-e, Bx-d))
    return {'theta3': t3, 'theta4': t4}

# Test Grashof condition for multiple configurations
configs = [
    {'a': 40, 'b': 120, 'c': 80, 'd': 100, 'e': 0, 'name': 'Crank-rocker'},
    {'a': 80, 'b': 120, 'c': 80, 'd': 100, 'e': 0, 'name': 'Non-Grashof'},
    {'a': 60, 'b': 120, 'c': 80, 'd': 100, 'e': 0, 'name': 'Change-point'},
    {'a': 60, 'b': 80, 'c': 70, 'd': 40, 'e': 0, 'name': 'Double-crank'},
]

print("\nGrashof Classification:")
for cfg in configs:
    D = np.sqrt(cfg['d']**2 + cfg['e']**2)
    links = sorted([cfg['a'], cfg['b'], cfg['c'], D])
    s, l = links[0], links[3]
    grashof = 'Grashof' if s + l <= links[1] + links[2] else 'Non-Grashof'
    if abs(s + l - links[1] - links[2]) < 0.01: grashof = 'Change-point'
    valid = sum(1 for t in range(361) if solve_fourbar(cfg['a'],cfg['b'],cfg['c'],cfg['d'],cfg['e'],t))
    print(f"  {cfg['name']}: s+l={s+l:.0f}, p+q={links[1]+links[2]:.0f} -> {grashof}, valid={valid}/361")

# ============================================================
# STEP 3: If CSV loaded, plot simulator data vs analytical
# ============================================================
if sim_data:
    sim_angles = [float(r['Input Angle (deg)']) for r in sim_data]
    sim_theta3 = [float(r['Coupler Angle (deg)']) for r in sim_data]
    sim_theta4 = [float(r['Follower Angle (deg)']) for r in sim_data]

    # Compute analytical for baseline config
    a, b, c, d, e = 40, 120, 80, 100, 0
    anal_theta3, anal_theta4 = [], []
    for t in sim_angles:
        r = solve_fourbar(a, b, c, d, e, t)
        anal_theta3.append(r['theta3'] if r else np.nan)
        anal_theta4.append(r['theta4'] if r else np.nan)

    fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
    axes[0].plot(sim_angles, sim_theta3, 'b-', label='Simulator', linewidth=2)
    axes[0].plot(sim_angles, anal_theta3, 'r--', label='Analytical', linewidth=1)
    axes[0].set_ylabel('Coupler Angle (deg)')
    axes[0].legend()
    axes[0].set_title('Experiment 1: Simulator vs Analytical Comparison')

    axes[1].plot(sim_angles, sim_theta4, 'b-', label='Simulator', linewidth=2)
    axes[1].plot(sim_angles, anal_theta4, 'r--', label='Analytical', linewidth=1)
    axes[1].set_ylabel('Follower Angle (deg)')
    axes[1].set_xlabel('Input Angle (deg)')
    axes[1].legend()

    plt.tight_layout()
    plt.savefig('plots/experiment_1_comparison.png', dpi=150)
    plt.show()

    # Compute error
    errors = [abs(s - a) for s, a in zip(sim_theta4, anal_theta4) if not np.isnan(a)]
    print(f"\nFollower angle: max error = {max(errors):.4f} deg, mean error = {np.mean(errors):.4f} deg")

import numpy as np
import matplotlib.pyplot as plt
import csv
import os

# ============================================================
# STEP 1: Load CSV data from the simulator
# Download CSV from the simulator or record from charts
# ============================================================
csv_files = {
    'Baseline': 'data/exp2_baseline.csv',
    'Short follower': 'data/exp2_short_follower.csv',
    'Improved': 'data/exp2_improved.csv',
}
sim_datasets = {}
for label, path in csv_files.items():
    if os.path.exists(path):
        with open(path) as f:
            sim_datasets[label] = list(csv.DictReader(f))
        print(f"Loaded {len(sim_datasets[label])} rows from {path}")
    else:
        print(f"No CSV at {path} (will use analytical only for {label})")

# ============================================================
# STEP 2: Analytical computation (two-circle intersection)
# ============================================================
def solve_fourbar(a, b, c, d, e, theta2_deg):
    t2 = np.radians(theta2_deg)
    Ax, Ay = a*np.cos(t2), a*np.sin(t2)
    dx, dy = d - Ax, e - Ay
    dist = np.sqrt(dx**2 + dy**2)
    if dist > b+c or dist < abs(b-c) or dist < 1e-10:
        return None
    ca = (b**2 - c**2 + dist**2) / (2*dist)
    hsq = b**2 - ca**2
    if hsq < 0: return None
    h = np.sqrt(hsq)
    mx, my = Ax + ca*dx/dist, Ay + ca*dy/dist
    Bx, By = mx + h*dy/dist, my - h*dx/dist
    t3 = np.degrees(np.arctan2(By-Ay, Bx-Ax))
    t4 = np.degrees(np.arctan2(By-e, Bx-d))
    mu = abs(t3 - t4)
    if mu > 180: mu = 360 - mu
    if mu > 90: mu = 180 - mu
    return {'theta3': t3, 'theta4': t4, 'mu': mu}

configs = [
    {'a': 40, 'b': 120, 'c': 80, 'd': 100, 'e': 0, 'label': 'Baseline'},
    {'a': 40, 'b': 120, 'c': 40, 'd': 100, 'e': 0, 'label': 'Short follower'},
    {'a': 30, 'b': 100, 'c': 90, 'd': 100, 'e': 0, 'label': 'Improved'},
]

theta2 = np.arange(0, 361)
fig, ax = plt.subplots(figsize=(10, 5))

# ============================================================
# STEP 3: Plot simulator data vs analytical
# ============================================================
for cfg in configs:
    # Analytical
    mu_anal = []
    for t in theta2:
        r = solve_fourbar(cfg['a'], cfg['b'], cfg['c'], cfg['d'], cfg['e'], t)
        mu_anal.append(r['mu'] if r else np.nan)
    mu_anal = np.array(mu_anal)

    label = cfg['label']
    if label in sim_datasets:
        sim = sim_datasets[label]
        sim_angles = [float(row['Input Angle (deg)']) for row in sim]
        sim_mu = [float(row['Transmission Angle (deg)']) for row in sim]
        ax.plot(sim_angles, sim_mu, '-', label=f'{label} (simulator)', linewidth=2)
        ax.plot(theta2, mu_anal, '--', label=f'{label} (analytical)', linewidth=1, alpha=0.7)
        # Error
        common = min(len(sim_mu), len(mu_anal))
        errors = [abs(s - a) for s, a in zip(sim_mu[:common], mu_anal[:common]) if not np.isnan(a)]
        if errors:
            print(f"{label}: max error = {max(errors):.4f} deg, mean = {np.mean(errors):.4f} deg")
    else:
        ax.plot(theta2, mu_anal, '-', label=label, linewidth=2)

    valid = mu_anal[~np.isnan(mu_anal)]
    print(f"{label}: min mu = {np.nanmin(mu_anal):.1f}, max mu = {np.nanmax(mu_anal):.1f}")

ax.axhline(y=40, color='r', linestyle='--', alpha=0.5, label='40 deg limit')
ax.set_xlabel('Input Angle (degrees)')
ax.set_ylabel('Transmission Angle (degrees)')
ax.set_title('Experiment 2: Transmission Angle Comparison')
ax.legend(fontsize=8)
plt.tight_layout()
plt.savefig('plots/experiment_2_transmission.png', dpi=150)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import csv
import os

# ============================================================
# STEP 1: Load CSV data from the simulator for each preset
# Download CSV from the simulator for each preset configuration
# ============================================================
csv_files = {
    'Crank-Rocker': 'data/exp3_crank_rocker.csv',
    'Double-Crank': 'data/exp3_double_crank.csv',
    'Double-Rocker': 'data/exp3_double_rocker.csv',
    'Parallelogram': 'data/exp3_parallelogram.csv',
}
sim_datasets = {}
for name, path in csv_files.items():
    if os.path.exists(path):
        with open(path) as f:
            sim_datasets[name] = list(csv.DictReader(f))
        print(f"Loaded {len(sim_datasets[name])} rows for {name}")
    else:
        print(f"No CSV at {path} (analytical only for {name})")

# ============================================================
# STEP 2: Analytical computation (two-circle intersection)
# ============================================================
def solve_fourbar(a, b, c, d, e, theta2_deg):
    t2 = np.radians(theta2_deg)
    Ax, Ay = a*np.cos(t2), a*np.sin(t2)
    dx, dy = d - Ax, e - Ay
    dist = np.sqrt(dx**2 + dy**2)
    if dist > b+c or dist < abs(b-c) or dist < 1e-10:
        return None
    ca = (b**2 - c**2 + dist**2) / (2*dist)
    hsq = b**2 - ca**2
    if hsq < 0: return None
    h = np.sqrt(hsq)
    mx, my = Ax + ca*dx/dist, Ay + ca*dy/dist
    Bx, By = mx + h*dy/dist, my - h*dx/dist
    t3 = np.degrees(np.arctan2(By-Ay, Bx-Ax))
    t4 = np.degrees(np.arctan2(By-e, Bx-d))
    return {'theta3': t3, 'theta4': t4}

presets = [
    {'name': 'Crank-Rocker', 'a': 40, 'b': 120, 'c': 80, 'd': 100, 'e': 0},
    {'name': 'Double-Crank', 'a': 60, 'b': 80, 'c': 70, 'd': 40, 'e': 0},
    {'name': 'Double-Rocker', 'a': 80, 'b': 35, 'c': 60, 'd': 100, 'e': 0},
    {'name': 'Parallelogram', 'a': 60, 'b': 100, 'c': 60, 'd': 100, 'e': 0},
]

theta2 = np.arange(0, 361)

# ============================================================
# STEP 3: Plot simulator data vs analytical for each preset
# ============================================================
fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
colors = ['blue', 'red', 'green', 'orange']

for idx, p in enumerate(presets):
    # Analytical
    t4_anal = []
    for t in theta2:
        r = solve_fourbar(p['a'], p['b'], p['c'], p['d'], p['e'], t)
        t4_anal.append(r['theta4'] if r else np.nan)
    t4_anal = np.array(t4_anal)

    name = p['name']
    if name in sim_datasets:
        sim = sim_datasets[name]
        sim_angles = [float(row['Input Angle (deg)']) for row in sim]
        sim_t4 = [float(row['Follower Angle (deg)']) for row in sim]
        sim_omega4 = [float(row['Omega4 (rad/s)']) for row in sim]
        axes[0].plot(sim_angles, sim_t4, '-', color=colors[idx],
                     label=f'{name} (sim)', linewidth=2)
        axes[0].plot(theta2, t4_anal, '--', color=colors[idx],
                     label=f'{name} (anal)', linewidth=1, alpha=0.6)
        axes[1].plot(sim_angles, sim_omega4, '-', color=colors[idx],
                     label=f'{name}', linewidth=2)
    else:
        axes[0].plot(theta2, t4_anal, '-', color=colors[idx],
                     label=name, linewidth=2)

    valid = t4_anal[~np.isnan(t4_anal)]
    if len(valid) > 0:
        osc = np.nanmax(t4_anal) - np.nanmin(t4_anal)
        print(f"{name}: theta4 = {np.nanmin(t4_anal):.1f} to {np.nanmax(t4_anal):.1f} ({osc:.1f} deg)")

axes[0].set_ylabel('Follower Angle theta4 (deg)')
axes[0].set_title('Experiment 3: Four-Bar Presets Comparison')
axes[0].legend(fontsize=8)
axes[1].set_ylabel('Follower Omega4 (rad/s)')
axes[1].set_xlabel('Input Angle theta2 (deg)')
axes[1].legend(fontsize=8)

plt.tight_layout()
plt.savefig('plots/experiment_3_presets.png', dpi=150)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import csv
import os

# ============================================================
# STEP 1: Load CSV data from the simulator
# Run the simulator in Open circuit, download CSV
# Then switch to Crossed circuit, download CSV
# ============================================================
csv_files = {
    'Open': 'data/exp4_open.csv',
    'Crossed': 'data/exp4_crossed.csv',
}
sim_datasets = {}
for label, path in csv_files.items():
    if os.path.exists(path):
        with open(path) as f:
            sim_datasets[label] = list(csv.DictReader(f))
        print(f"Loaded {len(sim_datasets[label])} rows for {label}")
    else:
        print(f"No CSV at {path} (analytical only for {label})")

# ============================================================
# STEP 2: Analytical computation (two-circle intersection)
# Open circuit: B above the line A-O4
# Crossed circuit: B below the line A-O4
# ============================================================
def solve_fourbar_circuit(a, b, c, d, e, theta2_deg, circuit='open'):
    t2 = np.radians(theta2_deg)
    Ax, Ay = a*np.cos(t2), a*np.sin(t2)
    dx, dy = d - Ax, e - Ay
    dist = np.sqrt(dx**2 + dy**2)
    if dist > b+c or dist < abs(b-c) or dist < 1e-10:
        return None
    ca = (b**2 - c**2 + dist**2) / (2*dist)
    hsq = b**2 - ca**2
    if hsq < 0: return None
    h = np.sqrt(hsq)
    mx, my = Ax + ca*dx/dist, Ay + ca*dy/dist
    sign = 1 if circuit == 'open' else -1
    Bx = mx + sign * h * dy / dist
    By = my - sign * h * dx / dist
    t4 = np.degrees(np.arctan2(By - e, Bx - d))
    t3 = np.degrees(np.arctan2(By - Ay, Bx - Ax))
    mu = abs(t3 - t4)
    if mu > 180: mu = 360 - mu
    if mu > 90: mu = 180 - mu
    return {'theta4': t4, 'mu': mu}

a, b, c, d, e = 40, 120, 80, 100, 0
theta2 = np.arange(0, 361)

# ============================================================
# STEP 3: Plot simulator data vs analytical
# ============================================================
fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
circuit_colors = {'Open': 'blue', 'Crossed': 'red'}

for circuit in ['Open', 'Crossed']:
    # Analytical
    t4_anal, mu_anal = [], []
    for t in theta2:
        r = solve_fourbar_circuit(a, b, c, d, e, t, circuit.lower())
        t4_anal.append(r['theta4'] if r else np.nan)
        mu_anal.append(r['mu'] if r else np.nan)
    t4_anal = np.array(t4_anal)
    mu_anal = np.array(mu_anal)
    color = circuit_colors[circuit]

    if circuit in sim_datasets:
        sim = sim_datasets[circuit]
        sim_angles = [float(row['Input Angle (deg)']) for row in sim]
        sim_t4 = [float(row['Follower Angle (deg)']) for row in sim]
        sim_mu = [float(row['Transmission Angle (deg)']) for row in sim]
        axes[0].plot(sim_angles, sim_t4, '-', color=color,
                     label=f'{circuit} (simulator)', linewidth=2)
        axes[0].plot(theta2, t4_anal, '--', color=color,
                     label=f'{circuit} (analytical)', linewidth=1, alpha=0.6)
        axes[1].plot(sim_angles, sim_mu, '-', color=color, label=f'{circuit} (sim)', linewidth=2)
        axes[1].plot(theta2, mu_anal, '--', color=color, label=f'{circuit} (anal)', linewidth=1, alpha=0.6)
    else:
        axes[0].plot(theta2, t4_anal, '-', color=color, label=circuit, linewidth=2)
        axes[1].plot(theta2, mu_anal, '-', color=color, label=circuit, linewidth=2)

    valid = t4_anal[~np.isnan(t4_anal)]
    if len(valid) > 0:
        print(f"{circuit}: theta4 = {np.nanmin(t4_anal):.1f} to {np.nanmax(t4_anal):.1f}")
        print(f"{circuit}: min mu = {np.nanmin(mu_anal):.1f}, max mu = {np.nanmax(mu_anal):.1f}")

axes[0].set_ylabel('Follower Angle (deg)')
axes[0].set_title('Experiment 4: Open vs Crossed Circuit')
axes[0].legend(fontsize=8)
axes[1].set_ylabel('Transmission Angle (deg)')
axes[1].set_xlabel('Input Angle (deg)')
axes[1].axhline(y=40, color='gray', linestyle='--', alpha=0.5, label='40 deg limit')
axes[1].legend(fontsize=8)

plt.tight_layout()
plt.savefig('plots/experiment_4_circuits.png', dpi=150)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import csv
import os

# ============================================================
# STEP 1: Load CSV data from the simulator
# The CSV includes Ax, Ay, Bx, By columns for coupler midpoint
# ============================================================
csv_files = {
    'Crank-Rocker': 'data/exp5_crank_rocker.csv',
    'Long coupler': 'data/exp5_long_coupler.csv',
    'Parallelogram': 'data/exp5_parallelogram.csv',
}
sim_datasets = {}
for name, path in csv_files.items():
    if os.path.exists(path):
        with open(path) as f:
            sim_datasets[name] = list(csv.DictReader(f))
        print(f"Loaded {len(sim_datasets[name])} rows for {name}")
    else:
        print(f"No CSV at {path} (analytical only for {name})")

# ============================================================
# STEP 2: Analytical computation (two-circle intersection)
# ============================================================
def solve_fourbar(a, b, c, d, e, theta2_deg):
    t2 = np.radians(theta2_deg)
    Ax, Ay = a*np.cos(t2), a*np.sin(t2)
    dx, dy = d - Ax, e - Ay
    dist = np.sqrt(dx**2 + dy**2)
    if dist > b+c or dist < abs(b-c) or dist < 1e-10:
        return None
    ca_val = (b**2 - c**2 + dist**2) / (2*dist)
    hsq = b**2 - ca_val**2
    if hsq < 0: return None
    h = np.sqrt(hsq)
    mx, my = Ax + ca_val*dx/dist, Ay + ca_val*dy/dist
    Bx, By = mx + h*dy/dist, my - h*dx/dist
    return {'Ax': Ax, 'Ay': Ay, 'Bx': Bx, 'By': By}

configs = [
    {'a': 40, 'b': 120, 'c': 80, 'd': 100, 'e': 0, 'name': 'Crank-Rocker'},
    {'a': 30, 'b': 150, 'c': 100, 'd': 100, 'e': 0, 'name': 'Long coupler'},
    {'a': 60, 'b': 100, 'c': 60, 'd': 100, 'e': 0, 'name': 'Parallelogram'},
]

theta2 = np.arange(0, 361)

# ============================================================
# STEP 3: Plot simulator coupler curves vs analytical
# ============================================================
fig, axes = plt.subplots(1, 3, figsize=(14, 5))

for idx, cfg in enumerate(configs):
    name = cfg['name']

    # Analytical coupler midpoint
    mx_a, my_a = [], []
    for t in theta2:
        r = solve_fourbar(cfg['a'], cfg['b'], cfg['c'], cfg['d'], cfg['e'], t)
        if r is None: continue
        mx_a.append((r['Ax'] + r['Bx']) / 2)
        my_a.append((r['Ay'] + r['By']) / 2)

    if name in sim_datasets:
        sim = sim_datasets[name]
        # Compute coupler midpoint from simulator joint coordinates
        sim_mx = [(float(row['Ax (mm)']) + float(row['Bx (mm)'])) / 2 for row in sim]
        sim_my = [(float(row['Ay (mm)']) + float(row['By (mm)'])) / 2 for row in sim]
        axes[idx].plot(sim_mx, sim_my, 'b-', label='Simulator', linewidth=2)
        axes[idx].plot(mx_a, my_a, 'r--', label='Analytical', linewidth=1)
        axes[idx].legend(fontsize=8)
    else:
        axes[idx].plot(mx_a, my_a, 'purple', linewidth=2)

    axes[idx].set_aspect('equal')
    axes[idx].set_title(name)
    axes[idx].set_xlabel('X (mm)')
    axes[idx].set_ylabel('Y (mm)')
    axes[idx].grid(True, alpha=0.3)

plt.suptitle('Experiment 5: Coupler Midpoint Curves')
plt.tight_layout()
plt.savefig('plots/experiment_5_coupler.png', dpi=150)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import csv
import os

# ============================================================
# STEP 1: Load CSV data from the simulator
# Run the simulator for each parameter variation, download CSV
# Name files: exp6_a30.csv, exp6_a40.csv, exp6_b100.csv, etc.
# ============================================================
csv_dir = 'data'
sim_data_a, sim_data_b, sim_data_d = {}, {}, {}

for a_val in [30, 35, 40, 45, 50, 55, 60]:
    path = os.path.join(csv_dir, f'exp6_a{a_val}.csv')
    if os.path.exists(path):
        with open(path) as f:
            rows = list(csv.DictReader(f))
        t4 = [float(r['Follower Angle (deg)']) for r in rows if float(r['Follower Angle (deg)']) != 0]
        if t4:
            sim_data_a[a_val] = max(t4) - min(t4)

for b_val in [100, 110, 120, 130, 140]:
    path = os.path.join(csv_dir, f'exp6_b{b_val}.csv')
    if os.path.exists(path):
        with open(path) as f:
            rows = list(csv.DictReader(f))
        t4 = [float(r['Follower Angle (deg)']) for r in rows if float(r['Follower Angle (deg)']) != 0]
        if t4:
            sim_data_b[b_val] = max(t4) - min(t4)

for d_val in [80, 90, 100, 110, 120]:
    path = os.path.join(csv_dir, f'exp6_d{d_val}.csv')
    if os.path.exists(path):
        with open(path) as f:
            rows = list(csv.DictReader(f))
        t4 = [float(r['Follower Angle (deg)']) for r in rows if float(r['Follower Angle (deg)']) != 0]
        if t4:
            sim_data_d[d_val] = max(t4) - min(t4)

loaded = len(sim_data_a) + len(sim_data_b) + len(sim_data_d)
print(f"Loaded simulator data for {loaded} configurations")
if loaded == 0:
    print("No CSV files found. Using analytical data only.")
    print("Download CSVs from the simulator for each parameter variation.")

# ============================================================
# STEP 2: Analytical computation (two-circle intersection)
# ============================================================
def get_oscillation(a, b, c, d, e=0):
    t4_arr = []
    for t in range(361):
        t2 = np.radians(t)
        Ax, Ay = a*np.cos(t2), a*np.sin(t2)
        dx, dy = d - Ax, e - Ay
        dist = np.sqrt(dx**2 + dy**2)
        if dist > b+c or dist < abs(b-c) or dist < 1e-10:
            continue
        ca = (b**2 - c**2 + dist**2) / (2*dist)
        hsq = b**2 - ca**2
        if hsq < 0: continue
        h = np.sqrt(hsq)
        mx, my = Ax + ca*dx/dist, Ay + ca*dy/dist
        Bx, By = mx + h*dy/dist, my - h*dx/dist
        t4_arr.append(np.degrees(np.arctan2(By-e, Bx-d)))
    if len(t4_arr) < 10: return None
    return max(t4_arr) - min(t4_arr)

# ============================================================
# STEP 3: Plot simulator data vs analytical
# ============================================================
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

# Vary a
a_vals = np.arange(25, 65, 5)
osc_a = [get_oscillation(a, 120, 80, 100) for a in a_vals]
axes[0].plot(a_vals, osc_a, 'b-', label='Analytical', linewidth=2)
if sim_data_a:
    sa_x = sorted(sim_data_a.keys())
    axes[0].plot(sa_x, [sim_data_a[k] for k in sa_x], 'ro', label='Simulator', markersize=8)
    axes[0].legend(fontsize=8)
axes[0].set_xlabel('Crank Length a (mm)')
axes[0].set_ylabel('Follower Oscillation (deg)')
axes[0].set_title('Sensitivity to a')

# Vary b
b_vals = np.arange(100, 145, 5)
osc_b = [get_oscillation(40, b, 80, 100) for b in b_vals]
axes[1].plot(b_vals, osc_b, 'r-', label='Analytical', linewidth=2)
if sim_data_b:
    sb_x = sorted(sim_data_b.keys())
    axes[1].plot(sb_x, [sim_data_b[k] for k in sb_x], 'bo', label='Simulator', markersize=8)
    axes[1].legend(fontsize=8)
axes[1].set_xlabel('Coupler Length b (mm)')
axes[1].set_title('Sensitivity to b')

# Vary d
d_vals = np.arange(75, 125, 5)
osc_d = [get_oscillation(40, 120, 80, d) for d in d_vals]
axes[2].plot(d_vals, osc_d, 'g-', label='Analytical', linewidth=2)
if sim_data_d:
    sd_x = sorted(sim_data_d.keys())
    axes[2].plot(sd_x, [sim_data_d[k] for k in sd_x], 'mo', label='Simulator', markersize=8)
    axes[2].legend(fontsize=8)
axes[2].set_xlabel('Ground Length d (mm)')
axes[2].set_title('Sensitivity to d')

plt.suptitle('Experiment 6: Parametric Sensitivity of Follower Oscillation')
plt.tight_layout()
plt.savefig('plots/experiment_6_sensitivity.png', dpi=150)
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import csv
import os

# ============================================================
# STEP 1: Load CSV data from the simulator
# Run at 60 RPM and 120 RPM, download CSV for each
# ============================================================
csv_files = {
    '60 RPM': 'data/exp7_60rpm.csv',
    '120 RPM': 'data/exp7_120rpm.csv',
}
sim_datasets = {}
for label, path in csv_files.items():
    if os.path.exists(path):
        with open(path) as f:
            sim_datasets[label] = list(csv.DictReader(f))
        print(f"Loaded {len(sim_datasets[label])} rows for {label}")
    else:
        print(f"No CSV at {path} (analytical only for {label})")

# ============================================================
# STEP 2: Analytical computation (two-circle intersection + velocity)
# ============================================================
def solve_fourbar_vel(a, b, c, d, e, theta2_deg, omega2):
    t2 = np.radians(theta2_deg)
    Ax, Ay = a*np.cos(t2), a*np.sin(t2)
    dx, dy = d - Ax, e - Ay
    dist = np.sqrt(dx**2 + dy**2)
    if dist > b+c or dist < abs(b-c) or dist < 1e-10:
        return None
    ca = (b**2 - c**2 + dist**2) / (2*dist)
    hsq = b**2 - ca**2
    if hsq < 0: return None
    h = np.sqrt(hsq)
    mx, my = Ax + ca*dx/dist, Ay + ca*dy/dist
    Bx, By = mx + h*dy/dist, my - h*dx/dist
    t3 = np.arctan2(By-Ay, Bx-Ax)
    t4 = np.arctan2(By-e, Bx-d)
    sin34 = np.sin(t3 - t4)
    if abs(sin34) < 1e-10: return None
    w3 = a*omega2*np.sin(t4 - t2) / (b*sin34)
    w4 = a*omega2*np.sin(t2 - t3) / (c*np.sin(t4 - t3))
    return {'omega3': w3, 'omega4': w4}

a, b, c, d, e_off = 40, 120, 80, 100, 0
theta2 = np.arange(0, 361)

# ============================================================
# STEP 3: Plot simulator data vs analytical
# ============================================================
fig, axes = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
rpm_colors = {'60 RPM': 'blue', '120 RPM': 'red'}

for rpm, color, label in [(60, 'blue', '60 RPM'), (120, 'red', '120 RPM')]:
    omega2 = rpm * 2 * np.pi / 60
    # Analytical
    w4_anal, w3_anal = [], []
    for t in theta2:
        r = solve_fourbar_vel(a, b, c, d, e_off, t, omega2)
        w4_anal.append(r['omega4'] if r else np.nan)
        w3_anal.append(r['omega3'] if r else np.nan)
    w4_anal = np.array(w4_anal)
    w3_anal = np.array(w3_anal)

    if label in sim_datasets:
        sim = sim_datasets[label]
        sim_angles = [float(row['Input Angle (deg)']) for row in sim]
        sim_w4 = [float(row['Omega4 (rad/s)']) for row in sim]
        sim_w3 = [float(row['Omega3 (rad/s)']) for row in sim]
        axes[0].plot(sim_angles, sim_w4, '-', color=color,
                     label=f'{label} (simulator)', linewidth=2)
        axes[0].plot(theta2, w4_anal, '--', color=color,
                     label=f'{label} (analytical)', linewidth=1, alpha=0.6)
        axes[1].plot(sim_angles, sim_w3, '-', color=color,
                     label=f'{label} (simulator)', linewidth=2)
        axes[1].plot(theta2, w3_anal, '--', color=color,
                     label=f'{label} (analytical)', linewidth=1, alpha=0.6)
        # Error report
        common = min(len(sim_w4), len(w4_anal))
        errs = [abs(s - a) for s, a in zip(sim_w4[:common], w4_anal[:common]) if not np.isnan(a)]
        if errs:
            print(f"{label}: omega4 max error = {max(errs):.4f} rad/s")
    else:
        axes[0].plot(theta2, w4_anal, '-', color=color, label=label, linewidth=2)
        axes[1].plot(theta2, w3_anal, '-', color=color, label=label, linewidth=2)

    print(f"{label}: peak |omega4| = {np.nanmax(np.abs(w4_anal)):.3f} rad/s")

axes[0].set_ylabel('Follower omega4 (rad/s)')
axes[0].set_title('Experiment 7: Angular Velocity at Different Speeds')
axes[0].legend(fontsize=8)
axes[0].axhline(y=0, color='k', linewidth=0.5)
axes[1].set_ylabel('Coupler omega3 (rad/s)')
axes[1].set_xlabel('Input Angle (degrees)')
axes[1].legend(fontsize=8)
axes[1].axhline(y=0, color='k', linewidth=0.5)

plt.tight_layout()
plt.savefig('plots/experiment_7_acceleration.png', dpi=150)
plt.show()

print(f"\nSpeed ratio: 2x -> omega4 ratio should be 2x (linear)")
print(f"Acceleration ratio should be 4x (quadratic)")