"""Turns a recorded run into horsepower, acceleration and timing numbers. Physics: the force pushing the bike forward at the rear tire is F = m_eff * a + m * g * (Crr + slope) + 1/2 * rho * CdA * v^2 (speeding up, + tire rolling resistance, + hill, + air drag). Power = F * v. Divide by 745.7 for horsepower at the wheel. Engine horsepower is the wheel figure divided by the drivetrain efficiency (clutch/chain or torque converter). Speed comes from the GPS (Doppler, ~10 per second) or from a wheel RPM sensor if one is fitted. Acceleration comes from an accelerometer if one is fitted, otherwise from how fast the smoothed speed changes. """ import math import numpy as np G = 9.81 RHO_AIR = 1.20 # kg/m^3, near sea level at room temperature LB_TO_KG = 0.45359237 MPS_TO_MPH = 2.236936 W_PER_HP = 745.7 M_TO_FT = 3.28084 # ------------------------------------------------------------ helpers --- def savgol(values, window, deriv=0, dt=1.0): """Savitzky-Golay filter (quadratic): smooths, or gives a smooth derivative.""" values = np.asarray(values, dtype=float) n = len(values) window = int(window) | 1 if n < 3: return np.zeros(n) if deriv else values.copy() if window > n: window = n if n % 2 else n - 1 if window < 3: return np.gradient(values, dt) if deriv else values.copy() half = window // 2 x = np.arange(-half, half + 1, dtype=float) basis = np.vander(x, 3, increasing=True) # 1, x, x^2 coeffs = np.linalg.pinv(basis)[deriv] # row 0 = value, row 1 = slope padded = np.concatenate([np.full(half, values[0]), values, np.full(half, values[-1])]) out = np.convolve(padded, coeffs[::-1], mode="valid") return out / dt if deriv else out def resample(times, values, grid): times = np.asarray(times, dtype=float) values = np.asarray(values, dtype=float) if len(times) < 2: return np.full(len(grid), values[0] if len(values) else 0.0) order = np.argsort(times) return np.interp(grid, times[order], values[order]) def crossing_time(t, v, level): """First time v reaches `level`, interpolated between samples, or None.""" above = np.nonzero(v >= level)[0] if not len(above): return None i = above[0] if i == 0: return float(t[0]) v0, v1 = v[i - 1], v[i] frac = (level - v0) / (v1 - v0) if v1 != v0 else 0.0 return float(t[i - 1] + frac * (t[i] - t[i - 1])) def efficiency(settings): if settings.get("drivetrain") == "torque_converter": return float(settings["efficiency_torque_converter"]) return float(settings["efficiency_clutch"]) def r(x, digits=2): return None if x is None or not math.isfinite(x) else round(float(x), digits) # ----------------------------------------------------------- analysis --- def analyze(run, settings): channels = run.get("channels", {}) warnings = [] mass = (float(settings["bike_weight_lb"]) + float(settings["rider_weight_lb"])) * LB_TO_KG mass_eff = mass * (1 + float(settings["rotating_pct"]) / 100.0) if not settings.get("weights_set"): warnings.append("Using example weights - enter your real weights in Settings.") # ---- speed (m/s) on an even time grid ---- wheel = channels.get("wheel_rpm") use_wheel = settings.get("speed_source") == "wheel" or ( settings.get("speed_source") == "auto" and wheel and max(wheel.get("v") or [0]) > 0) gps = channels.get("gps") or {} if use_wheel and wheel and len(wheel.get("t", [])) >= 5: source = "wheel sensor" t_raw = np.asarray(wheel["t"], dtype=float) circumference = math.pi * float(settings["tire_diameter_in"]) * 0.0254 v_raw = np.asarray(wheel["v"], dtype=float) / 60.0 * circumference dt = 0.05 else: source = "GPS" ok = np.asarray(gps.get("ok", []), dtype=bool) t_raw = np.asarray(gps.get("t", []), dtype=float) v_raw = np.asarray(gps.get("speed_mps", []), dtype=float) if len(ok) == len(t_raw) and len(t_raw): if ok.mean() < 0.9: warnings.append("GPS fix was weak for part of this run - park in open sky first.") t_raw, v_raw = t_raw[ok], v_raw[ok] dt = float(np.median(np.diff(t_raw))) if len(t_raw) > 2 else 0.1 dt = min(max(dt, 0.02), 0.5) if dt > 0.15: warnings.append("GPS sent only %.0f updates a second (10 is best)." % (1 / dt)) sats = gps.get("sats") or [] if sats and np.mean(sats) < 8: warnings.append("Only %.0f satellites on average - speed is less exact." % np.mean(sats)) if len(t_raw) < 8: return {"ok": False, "error": "Not enough speed readings in this run.", "warnings": warnings} grid = np.arange(t_raw.min(), t_raw.max() + 1e-9, dt) window = int(settings.get("smoothing", 7)) if source == "wheel sensor": window = window * 2 + 1 # 20 readings/s instead of 10: same smoothing in seconds v = np.clip(savgol(resample(t_raw, v_raw, grid), window), 0, None) a = savgol(resample(t_raw, v_raw, grid), window, deriv=1, dt=dt) accel_source = "speed change (" + source + ")" accel = channels.get("accel_g") if accel and len(accel.get("t", [])) >= 10: a = savgol(resample(accel["t"], accel["v"], grid), window) * G accel_source = "accelerometer" # ---- launch and the "pull" (from launch to top speed) ---- t_launch = crossing_time(grid, v, 1.0 / MPS_TO_MPH) # 1 mph if t_launch is None: return {"ok": False, "error": "The bike never got moving in this run.", "warnings": warnings} i_launch = int(np.searchsorted(grid, t_launch)) i_top = i_launch + int(np.argmax(v[i_launch:])) pull = np.arange(len(grid)) pull = (pull >= i_launch) & (pull <= i_top) & (v > 2.0 / MPS_TO_MPH) # ---- forces, power ---- slope = float(settings["slope_pct"]) / 100.0 angle = math.atan(slope) f_roll = mass * G * float(settings["crr"]) * math.cos(angle) f_hill = mass * G * math.sin(angle) f_air = 0.5 * RHO_AIR * float(settings["cda_m2"]) * v ** 2 force = mass_eff * a + f_roll + f_hill + f_air hp_wheel = force * v / W_PER_HP hp_wheel_s = savgol(hp_wheel, 5) eff = efficiency(settings) hp_engine_s = hp_wheel_s / eff t = grid - t_launch mph = v * MPS_TO_MPH distance_ft = np.concatenate([[0.0], np.cumsum((v[1:] + v[:-1]) / 2 * np.diff(grid))]) * M_TO_FT distance_ft -= np.interp(t_launch, grid, distance_ft) summary = {"speed_source": source, "accel_source": accel_source, "mass_lb": r(mass / LB_TO_KG, 0), "efficiency": eff} if pull.sum() >= 3: hp_pull = np.where(pull, hp_wheel_s, -np.inf) i_peak = int(np.argmax(hp_pull)) summary.update({ "peak_wheel_hp": r(hp_wheel_s[i_peak]), "peak_engine_hp": r(hp_engine_s[i_peak]), "peak_hp_at_mph": r(mph[i_peak], 1), "max_accel_g": r(np.max(a[pull]) / G, 2), }) else: warnings.append("The run was too short to measure power.") summary.update({ "top_mph": r(mph[i_top], 1), "time_to_top_s": r(grid[i_top] - t_launch, 2), "zero_to_10_s": r(_since(crossing_time(grid, v, 10 / MPS_TO_MPH), t_launch)), "zero_to_20_s": r(_since(crossing_time(grid, v, 20 / MPS_TO_MPH), t_launch)), "zero_to_30_s": r(_since(crossing_time(grid, v, 30 / MPS_TO_MPH), t_launch)), "sixty_ft_s": r(_since(crossing_time(grid, distance_ft, 60.0), 0.0) if distance_ft[-1] >= 60 else None), "distance_ft": r(float(np.max(distance_ft)), 0), }) # ---- engine RPM: the classic dyno chart (power and torque by RPM) ---- by_rpm = None engine = channels.get("engine_rpm") if engine and len(engine.get("t", [])) >= 10 and max(engine["v"]) > 500: rpm = savgol(resample(engine["t"], engine["v"], grid), window) ok_rpm = pull & (rpm > 1000) if ok_rpm.sum() >= 3: torque = hp_engine_s * 5252.0 / np.maximum(rpm, 1) by_rpm = {"rpm": _list(rpm[ok_rpm], 0), "hp": _list(hp_engine_s[ok_rpm]), "torque_lbft": _list(torque[ok_rpm])} i = int(np.argmax(np.where(ok_rpm, hp_engine_s, -np.inf))) summary["peak_hp_at_rpm"] = r(rpm[i], 0) summary["peak_torque_lbft"] = r(float(np.max(torque[ok_rpm]))) keep = (t >= -1.0) series = { "t": _list(t[keep]), "mph": _list(mph[keep], 1), "accel_g": _list(a[keep] / G, 3), "hp_wheel": _list(np.where(v[keep] > 1 / MPS_TO_MPH, hp_wheel_s[keep], 0)), "distance_ft": _list(distance_ft[keep], 1), } for name in ("wheel_rpm", "engine_rpm"): ch = channels.get(name) if ch and len(ch.get("t", [])) >= 2: series[name] = _list(resample(ch["t"], ch["v"], grid)[keep], 0) curve = {"mph": _list(mph[pull], 1), "hp_wheel": _list(hp_wheel_s[pull]), "hp_engine": _list(hp_engine_s[pull])} return {"ok": True, "summary": summary, "series": series, "curve_by_speed": curve, "curve_by_rpm": by_rpm, "warnings": warnings} def _since(t, t0): return None if t is None else t - t0 def _list(arr, digits=2): return [round(float(x), digits) for x in arr] # --------------------------------------------------------- simulation --- def simulate_run(engine_hp=6.5, seconds=9.0, settings=None, top_mph=36.0, dt_out=0.1): """A believable mini bike pull, used by the practice GPS and the self-test. The engine makes `engine_hp` at the crank; the clutch/chain passes on the usual share; traction limits the launch; power fades near the governed top speed. Then it brakes to a stop. """ from .config import DEFAULT_SETTINGS s = dict(DEFAULT_SETTINGS) s.update(settings or {}) mass = (s["bike_weight_lb"] + s["rider_weight_lb"]) * LB_TO_KG mass_eff = mass * (1 + s["rotating_pct"] / 100.0) p_wheel = engine_hp * W_PER_HP * efficiency(s) f_traction = 0.45 * mass * G v_top = top_mph / MPS_TO_MPH angle = math.atan(s["slope_pct"] / 100.0) v, t, h = 0.0, 0.0, 0.002 out_t, out_v, next_out = [], [], 0.0 braking = False while True: if not braking and (t >= seconds or v >= v_top * 0.995): braking = True if braking: drive = -0.45 * mass * G else: power = p_wheel * max(0.0, 1 - (v / v_top) ** 8) drive = min(f_traction, power / max(v, 0.1)) resist = mass * G * (s["crr"] * math.cos(angle) + math.sin(angle)) + 0.5 * RHO_AIR * s["cda_m2"] * v * v v = max(0.0, v + (drive - resist) / mass_eff * h) t += h if t >= next_out: out_t.append(round(t, 3)) out_v.append(v) next_out += dt_out if braking and v <= 0.0: break return {"t": out_t, "speed_mps": out_v, "wheel_hp": p_wheel / W_PER_HP}