Public
Prueba del motor (code)
Code written by the AI agents of “Prueba del motor”. Each agent works on its own branch; main follows the lead agent.
Commitsmain.pyModified+61−45
RK4 projectile with quadratic drag, ground-crossing interpolation, no-drag analytic check
Lasques agent 1 (lead, iteration 1) · Oct 24, 2026 · a4e849d6a5f176997b6568f46823059f064890b2
1 file changed
+61−45| @@ -1,4 +1,3 @@ | ||
| 1 | """RK4 simulation of 2D projectile with quadratic drag a = -g y - (k/m)|v|v.""" | |
| 2 | 1 | import numpy as np |
| 3 | 2 | import lasques |
| 4 | 3 | import matplotlib |
| @@ -8,62 +7,79 @@ import matplotlib.pyplot as plt | ||
| 8 | 7 | p = lasques.params() |
| 9 | 8 | v0 = float(p.get("v0", 30.0)) |
| 10 | 9 | angle_deg = float(p.get("angle_deg", 45.0)) |
| 10 | g = float(p.get("g", 9.81)) | |
| 11 | 11 | km = float(p.get("km", 0.0)) |
| 12 | dt = float(p.get("dt", 1e-3)) | |
| 13 | g = 9.81 | |
| 12 | dt = float(p.get("dt", 0.001)) | |
| 14 | 13 | |
| 15 | th = np.radians(angle_deg) | |
| 16 | state = np.array([0.0, 0.0, v0*np.cos(th), v0*np.sin(th)]) | |
| 14 | theta = np.radians(angle_deg) | |
| 17 | 15 | |
| 18 | def deriv(s): | |
| 19 | x, y, vx, vy = s | |
| 16 | def accel(state, km): | |
| 17 | x, y, vx, vy = state | |
| 20 | 18 | sp = np.hypot(vx, vy) |
| 21 | return np.array([vx, vy, -km*sp*vx, -g - km*sp*vy]) | |
| 19 | return np.array([vx, vy, -km * sp * vx, -g - km * sp * vy]) | |
| 22 | 20 | |
| 21 | def rk4_step(state, dt, km): | |
| 22 | k1 = accel(state, km) | |
| 23 | k2 = accel(state + 0.5 * dt * k1, km) | |
| 24 | k3 = accel(state + 0.5 * dt * k2, km) | |
| 25 | k4 = accel(state + dt * k3, km) | |
| 26 | return state + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4) | |
| 27 | ||
| 28 | state = np.array([0.0, 0.0, v0 * np.cos(theta), v0 * np.sin(theta)]) | |
| 23 | 29 | traj = [state.copy()] |
| 24 | 30 | t = 0.0 |
| 25 | max_steps = 2_000_000 | |
| 31 | max_steps = 10_000_000 | |
| 26 | 32 | for i in range(max_steps): |
| 27 | k1 = deriv(state) | |
| 28 | k2 = deriv(state + 0.5*dt*k1) | |
| 29 | k3 = deriv(state + 0.5*dt*k2) | |
| 30 | k4 = deriv(state + dt*k3) | |
| 31 | new = state + dt/6*(k1 + 2*k2 + 2*k3 + k4) | |
| 32 | # detect ground crossing between state and new | |
| 33 | if new[1] < 0: | |
| 34 | # linear interpolation for range | |
| 35 | f = state[1] / (state[1] - new[1]) | |
| 36 | x_land = state[0] + f*(new[0]-state[0]) | |
| 37 | t_land = t + f*dt | |
| 38 | state = new | |
| 39 | traj.append(new) | |
| 40 | break | |
| 41 | state = new | |
| 42 | traj.append(new) | |
| 33 | prev = state | |
| 34 | state = rk4_step(state, dt, km) | |
| 43 | 35 | t += dt |
| 36 | traj.append(state.copy()) | |
| 37 | if state[1] < 0.0: # crossed ground: linear interpolation of last step | |
| 38 | f = prev[1] / (prev[1] - state[1]) | |
| 39 | range_m = prev[0] + f * (state[0] - prev[0]) | |
| 40 | break | |
| 41 | else: | |
| 42 | range_m = state[0] | |
| 44 | 43 | |
| 45 | 44 | traj = np.array(traj) |
| 46 | range_m = x_land | |
| 47 | apex = traj[:,1].max() | |
| 48 | R_analytic_no_drag = v0**2*np.sin(2*th)/g | |
| 49 | 45 | |
| 50 | lasques.set_metrics({"range_m": range_m, "flight_time_s": t_land, | |
| 51 | "apex_m": apex, "km": km, "dt": dt, | |
| 52 | "R_analytic_nodrag": R_analytic_no_drag}) | |
| 46 | # analytic no-drag reference | |
| 47 | R_analytic = v0**2 * np.sin(2 * theta) / g | |
| 48 | abs_drop = R_analytic - range_m if km == 0.0 else None | |
| 49 | ||
| 50 | lasques.set_metrics({ | |
| 51 | "range_m": range_m, | |
| 52 | "range_analytic_m": R_analytic, | |
| 53 | "km": km, | |
| 54 | "dt": dt, | |
| 55 | "angle_deg": angle_deg, | |
| 56 | "v0": v0, | |
| 57 | }) | |
| 53 | 58 | |
| 54 | if km == 0: | |
| 55 | rel_err = abs(range_m - R_analytic_no_drag)/R_analytic_no_drag | |
| 56 | lasques.set_metrics({"nodrag_rel_err": rel_err}) | |
| 57 | lasques.log_metric("nodrag_rel_err", rel_err) | |
| 58 | lasques.log_metric("range", range_m) | |
| 59 | lasques.save_table(traj[::max(1, len(traj)//2000), :], "trajectory") | |
| 59 | if km == 0.0: | |
| 60 | err_pct = abs(range_m - R_analytic) / R_analytic * 100 | |
| 61 | lasques.set_metrics({"no_drag_error_pct": err_pct}) | |
| 62 | lasques.summary( | |
| 63 | f"No-drag validation: simulated range {range_m:.4f} m vs analytic " | |
| 64 | f"{R_analytic:.4f} m (error {err_pct:.4e}%), dt={dt}, v0={v0}, angle={angle_deg}°" | |
| 65 | ) | |
| 66 | else: | |
| 67 | lasques.summary( | |
| 68 | f"k/m={km}: range {range_m:.4f} m at {angle_deg}°, v0={v0} m/s, dt={dt}" | |
| 69 | ) | |
| 60 | 70 | |
| 61 | fig, ax = plt.subplots(figsize=(7,4)) | |
| 62 | ax.plot(traj[:,0], traj[:,1], label=f"k/m={km}") | |
| 63 | ax.set_xlabel("x (m)"); ax.set_ylabel("y (m)") | |
| 64 | ax.set_title(f"Projectile, v0={v0} m/s, angle={angle_deg} deg, k/m={km}, dt={dt}") | |
| 65 | ax.legend(); ax.grid(True) | |
| 66 | lasques.save_figure(fig, f"traj_km{km}_dt{dt}") | |
| 71 | # trajectory figure | |
| 72 | fig, ax = plt.subplots(figsize=(7, 4)) | |
| 73 | ax.plot(traj[:, 0], traj[:, 1], lw=1.5, label=f"k/m={km}") | |
| 74 | ax.set_xlabel("x (m)") | |
| 75 | ax.set_ylabel("y (m)") | |
| 76 | ax.set_title(f"Trajectory, v0={v0} m/s, θ={angle_deg}°, dt={dt}") | |
| 77 | ax.legend() | |
| 78 | ax.grid(alpha=0.3) | |
| 79 | lasques.save_figure(fig, "trajectory") | |
| 67 | 80 | |
| 68 | lasques.summary(f"k/m={km}, dt={dt}: range={range_m:.4f} m, apex={apex:.3f} m, flight time={t_land:.3f} s") | |
| 69 | print(f"km={km} dt={dt} range={range_m:.5f} apex={apex:.4f} t={t_land:.4f}") | |
| 81 | lasques.save_table( | |
| 82 | {"t": np.arange(len(traj)) * dt, "x": traj[:, 0], "y": traj[:, 1], | |
| 83 | "vx": traj[:, 2], "vy": traj[:, 3]}, | |
| 84 | "trajectory", | |
| 85 | ) | |