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.

Commits

Single-run RK4 drag simulation with ground-crossing interpolation, per method 6737b2d4

Lasques agent 1 (lead, iteration 1) · Oct 24, 2026 · 49d58edebfeb79c9da81bf43e97369441ca93e4f

1 file changed

+57−134
main.pyModified+57−134
@@ -1,144 +1,67 @@
1"""Full k/m sweep + convergence + angle sensitivity in one run.
2Reuses the RK4 protocol from knowledge item 6737b2d4 (a = -g y_hat - (k/m)|v|v,
3RK4, linear interpolation of ground crossing)."""
1"""RK4 projectile with quadratic drag: k/m sweep, convergence check, analytic validation.
2Reuses method 6737b2d4-0136-463f-bc00-1e0e9a26bb10. One run = one (km, dt) combo."""
43import math
54import numpy as np
6import matplotlib
7matplotlib.use("Agg")
8import matplotlib.pyplot as plt
95import lasques
106
117p = lasques.params()
12G = 9.81
13V0 = 30.0
14DT = float(p.get("dt", 1e-3))
8g = float(p.get("g", 9.81))
9v0 = float(p.get("v0", 30.0))
10angle_deg = float(p.get("angle_deg", 45.0))
11km = float(p.get("km", 0.0))
12dt = float(p.get("dt", 1e-3))
13
14th = math.radians(angle_deg)
15s = np.array([0.0, 0.0, v0 * math.cos(th), v0 * math.sin(th)])
16
17
18def deriv(s):
19 vx, vy = s[2], s[3]
20 sp = math.hypot(vx, vy)
21 return np.array([vx, vy, -km * sp * vx, -g - km * sp * vy])
22
23
24def rk4_step(s, dt):
25 k1 = deriv(s)
26 k2 = deriv(s + 0.5 * dt * k1)
27 k3 = deriv(s + 0.5 * dt * k2)
28 k4 = deriv(s + dt * k3)
29 return s + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
30
31
32t = 0.0
33traj = []
34while t < 60.0:
35 traj.append((t, s[0], s[1], s[2], s[3]))
36 s_new = rk4_step(s, dt)
37 if s[1] > 0.0 >= s_new[1]:
38 f = s[1] / (s[1] - s_new[1]) # linear interpolation to ground
39 s_hit = s + f * (s_new - s)
40 t_hit = t + f * dt
41 break
42 s = s_new
43 t += dt
44else:
45 raise RuntimeError("no ground crossing within 60 s")
46
47R = float(s_hit[0])
48R_analytic = v0 ** 2 * math.sin(2 * th) / g
49err_vs_analytic = abs(R - R_analytic) / R_analytic if km == 0.0 else float("nan")
1550
16
17def deriv(s, km):
18 x, y, vx, vy = s
19 v = math.hypot(vx, vy)
20 return np.array([vx, vy, -km * v * vx, -G - km * v * vy])
21
22
23def simulate(km, angle_deg, dt, v0=V0):
24 th = math.radians(angle_deg)
25 s = np.array([0.0, 0.0, v0 * math.cos(th), v0 * math.sin(th)])
26 traj = [(0.0, 0.0)]
27 t = 0.0
28 while True:
29 k1 = deriv(s, km)
30 k2 = deriv(s + 0.5 * dt * k1, km)
31 k3 = deriv(s + 0.5 * dt * k2, km)
32 k4 = deriv(s + dt * k3, km)
33 s_new = s + (dt / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
34 t += dt
35 if s_new[1] < 0.0: # ground crossing: linear interpolation
36 f = s[1] / (s[1] - s_new[1])
37 s_hit = s + f * (s_new - s)
38 traj.append((t - dt + f * dt, s_hit[0]))
39 return s_hit[0], t - dt + f * dt, np.array(traj)
40 s = s_new
41 traj.append((t, s[0]))
42 if t > 100:
43 raise RuntimeError("no ground crossing")
44
45
46# --- 1. main sweep ---
47KM = [0.0, 0.0125, 0.025, 0.0375, 0.05]
48rows, trajs = [], {}
49for km in KM:
50 R, tf, traj = simulate(km, 45.0, DT)
51 trajs[km] = traj
52 rows.append({"km": km, "range_m": R, "flight_time_s": tf})
53
54R0 = rows[0]["range_m"]
55for r in rows:
56 r["drop_m"] = R0 - r["range_m"]
57 r["drop_pct"] = 100.0 * (R0 - r["range_m"]) / R0
58
59R_analytic = V0 ** 2 * math.sin(2 * math.radians(45.0)) / G
60nodrag_err_pct = 100 * abs(R0 - R_analytic) / R_analytic
61
62# --- 2. dt convergence ---
63conv = []
64for km in (0.0, 0.05):
65 prev = None
66 for dt in (4e-3, 2e-3, 1e-3, 5e-4, 2.5e-4):
67 R, _, _ = simulate(km, 45.0, dt)
68 rel = None if prev is None else abs(R - prev) / prev
69 conv.append({"km": km, "dt": dt, "range_m": R, "rel_change_vs_prev": rel})
70 prev = R
71conv_change_1e3_vs_5e4 = [
72 c for c in conv if c["dt"] in (1e-3, 5e-4)
73]
74max_conv_change = max(
75 abs(a["range_m"] - b["range_m"]) / b["range_m"]
76 for a in conv if a["dt"] == 1e-3 for b in conv if b["dt"] == 5e-4 and a["km"] == b["km"]
77)
78
79# --- 3. curve shape: fractional loss vs km, second differences ---
80frac = [r["drop_pct"] for r in rows]
81second_diffs = [frac[i] - 2 * frac[i - 1] + frac[i - 2] for i in range(2, len(frac))]
82
83# --- 4. optimal launch angle vs km ---
84opt_rows = []
85for km in (0.0, 0.0125, 0.025, 0.0375, 0.05):
86 best = max(
87 ((simulate(km, a, 1e-3)[0], a) for a in np.arange(20.0, 70.01, 0.5))
88 )
89 R45 = simulate(km, 45.0, 1e-3)[0]
90 opt_rows.append({"km": km, "opt_angle_deg": best[1], "opt_range_m": best[0],
91 "range_at_45": R45,
92 "gain_from_45_pct": 100 * (best[0] - R45) / R45})
93
94# --- outputs ---
9551lasques.set_metrics({
96 "nodrag_range_m": R0,
97 "nodrag_analytic_m": R_analytic,
98 "nodrag_err_pct": nodrag_err_pct,
99 "range_km0": rows[0]["range_m"],
100 "range_km0p0125": rows[1]["range_m"],
101 "range_km0p025": rows[2]["range_m"],
102 "range_km0p0375": rows[3]["range_m"],
103 "range_km0p05": rows[4]["range_m"],
104 "drop_pct_km0p0125": rows[1]["drop_pct"],
105 "drop_pct_km0p025": rows[2]["drop_pct"],
106 "drop_pct_km0p0375": rows[3]["drop_pct"],
107 "drop_pct_km0p05": rows[4]["drop_pct"],
108 "max_conv_rel_change_1e3_vs_5e4": max_conv_change,
109 "opt_angle_km0": opt_rows[0]["opt_angle_deg"],
110 "opt_angle_km0p05": opt_rows[4]["opt_angle_deg"],
52 "range_m": R,
53 "time_of_flight_s": float(t_hit),
54 "km": km,
55 "dt": dt,
56 "analytic_nodrag_range_m": R_analytic,
57 "nodrag_validation_relerr": err_vs_analytic,
11158})
112
113lasques.save_table(rows, "sweep_range_vs_km")
114lasques.save_table(conv, "dt_convergence")
115lasques.save_table(opt_rows, "optimal_angle_vs_km")
116
117fig, ax = plt.subplots(1, 2, figsize=(11, 4.2))
118for km, traj in trajs.items():
119 tr = np.array(traj)
120 ax[0].plot(tr[:, 1], tr[:, 0], label=f"k/m={km}")
121ax[0].set_xlabel("time (s)"); ax[0].set_ylabel("x (m)"); ax[0].legend(fontsize=7)
122ax[0].set_title("x(t) per drag coefficient")
123ax[1].plot([r["km"] for r in rows], [r["drop_pct"] for r in rows], "o-")
124ax[1].set_xlabel("k/m (1/m)"); ax[1].set_ylabel("range drop vs no drag (%)")
125ax[1].set_title("Range drop vs k/m (v0=30 m/s, 45°)")
126lasques.save_figure(fig, "drag_sweep")
127
59lasques.save_table(
60 [(r[0], r[1], r[2], r[3], r[4]) for r in traj[:: max(1, int(0.05 / dt))]],
61 "trajectory",
62)
12863lasques.summary(
129 "**Sweep (v0=30 m/s, 45°, RK4 dt=1e-3):**\n\n"
130 + "| k/m | Range (m) | Drop (m) | Drop % |\n|---|---|---|---|\n"
131 + "\n".join(
132 f"| {r['km']:.4f} | {r['range_m']:.3f} | {r['drop_m']:.3f} | {r['drop_pct']:.2f}% |"
133 for r in rows
134 )
135 + f"\n\nNo-drag vs analytic (91.77 m): error {nodrag_err_pct:.4f}%. "
136 + f"Max relative range change dt 1e-3 vs 5e-4: {max_conv_change:.2e}.\n\n"
137 + "Second differences of fractional loss vs k/m: "
138 + ", ".join(f"{d:.2f}" for d in second_diffs)
139 + " (all negative -> concave, saturating).\n\n"
140 + "Optimal launch angle: "
141 + "; ".join(f"k/m={o['km']:.4f}: {o['opt_angle_deg']:.1f}° (gain {o['gain_from_45_pct']:.2f}% vs 45°)"
142 for o in opt_rows)
64 f"km={km}, dt={dt}: range R = {R:.4f} m, t_flight = {t_hit:.4f} s. "
65 f"Analytic no-drag R = {R_analytic:.4f} m"
66 + (f" (rel err {err_vs_analytic:.2e})" if km == 0.0 else "")
14367)
144print("done")