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.

RK4 projectile with quadratic drag: single-run script and aggregate sweep script

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

2 files changed

+128−69

main.pyModified+64−69

@@ -1,85 +1,80 @@
	1	"""Aggregate run: sweep km values, convergence checks, plots, final table."""
1	2	import numpy as np
2	3	import lasques
3	4	import matplotlib
4	5	matplotlib.use("Agg")
5	6	import matplotlib.pyplot as plt
6	7
7		p = lasques.params()
8		v0 = float(p.get("v0", 30.0))
9		angle_deg = float(p.get("angle_deg", 45.0))
10		g = float(p.get("g", 9.81))
11		km = float(p.get("km", 0.0))
12		dt = float(p.get("dt", 0.001))
	8	g = 9.81
	9	v0 = 30.0
	10	angle_deg = 45.0
	11	km_values = [0.0, 0.0125, 0.025, 0.0375, 0.05]
	12	dt_values = [1e-2, 5e-3, 1e-3, 5e-4, 1e-4]
13	13
14		theta = np.radians(angle_deg)
	14	th = np.radians(angle_deg)
	15	R_analytic = v0*2np.sin(2*th)/g # no-drag analytic range
15	16
16		def accel(state, km):
17		x, y, vx, vy = state
18		sp = np.hypot(vx, vy)
19		return np.array([vx, vy, -km * sp * vx, -g - km * sp * vy])
	17	def simulate(km, dt):
	18	s = np.array([0.0, 0.0, v0np.cos(th), v0np.sin(th)])
	19	def deriv(s):
	20	x, y, vx, vy = s
	21	spd = np.hypot(vx, vy)
	22	return np.array([vx, vy, -kmspdvx, -g - kmspdvy])
	23	t = 0.0
	24	traj = [s.copy()]
	25	while s[1] >= 0.0:
	26	k1 = deriv(s); k2 = deriv(s+0.5dtk1); k3 = deriv(s+0.5dtk2); k4 = deriv(s+dt*k3)
	27	s_new = s + dt(k1+2k2+2*k3+k4)/6.0
	28	if s_new[1] < 0:
	29	frac = s[1]/(s[1]-s_new[1]); s = s + frac(s_new-s); t += fracdt
	30	traj.append(s); break
	31	s = s_new; t += dt; traj.append(s.copy())
	32	if t > 60: raise RuntimeError
	33	return np.array(traj), s[0], t
20	34
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)
	35	# Validation: no-drag across dt
	36	lasques.save_json({"R_analytic": R_analytic}, "analytic")
	37	for dt in dt_values:
	38	_, R, _ = simulate(0.0, dt)
	39	err = abs(R - R_analytic)/R_analytic*100
	40	lasques.set_metrics({f"nodrag_R_dt{dt}": R, f"nodrag_errpct_dt{dt}": err})
27	41
28		state = np.array([0.0, 0.0, v0 * np.cos(theta), v0 * np.sin(theta)])
29		traj = [state.copy()]
30		t = 0.0
31		max_steps = 10_000_000
32		for i in range(max_steps):
33		prev = state
34		state = rk4_step(state, dt, km)
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]
	42	# Main sweep at dt=1e-3 with convergence check at dt=1e-4
	43	rows = []
	44	trajs = {}
	45	for km in km_values:
	46	traj, R, tof = simulate(km, 1e-3)
	47	_, R_fine, _ = simulate(km, 1e-4)
	48	trajs[km] = traj
	49	apex = traj[:,1].max()
	50	rows.append({"km": km, "range_m": R, "tof_s": tof, "apex_m": apex,
	51	"range_dt1e-4": R_fine, "conv_err_pct": abs(R-R_fine)/R_fine*100})
43	52
44		traj = np.array(traj)
45
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
	53	R0 = rows[0]["range_m"]
	54	for r in rows:
	55	r["drop_m"] = R0 - r["range_m"]
	56	r["drop_pct"] = (R0 - r["range_m"])/R0*100
49	57
	58	import pandas as pd
	59	df = pd.DataFrame(rows)
	60	lasques.save_table(df, "results")
50	61	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,
	62	"R_nodrag": R0, "R_analytic": R_analytic,
	63	"validation_err_pct": abs(R0-R_analytic)/R_analytic*100,
	64	**{f"R_km{km}": r["range_m"] for km, r in zip(km_values, rows)},
	65	**{f"drop_pct_km{km}": r["drop_pct"] for km, r in zip(km_values, rows)},
57	66	})
58	67
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		)
70
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")
	68	fig, ax = plt.subplots(figsize=(7,4.5))
	69	for km in km_values:
	70	ax.plot(trajs[km][:,0], trajs[km][:,1], label=f"k/m={km}")
	71	ax.set_xlabel("x (m)"); ax.set_ylabel("y (m)"); ax.legend(); ax.grid(alpha=0.3)
	72	ax.set_title(f"Projectile trajectories, v0={v0} m/s, θ={angle_deg}°")
	73	lasques.save_figure(fig, "trajectories")
80	74
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		)
	75	md = [f"v0={v0} m/s, θ={angle_deg}°, RK4 dt=1e-3 s", "",
	76	"\| k/m \| range (m) \| drop (m) \| drop (%) \| conv err % \|",
	77	"\|---\|---\|---\|---\|---\|"]
	78	for r in rows:
	79	md.append(f"\| {r['km']} \| {r['range_m']:.3f} \| {r['drop_m']:.3f} \| {r['drop_pct']:.2f} \| {r['conv_err_pct']:.4f} \|")
	80	lasques.summary("\n".join(md))

projectile.pyAdded+64−0

@@ -0,0 +1,64 @@
	1	"""RK4 simulation of projectile motion with quadratic drag.
	2	State: (x, y, vx, vy). a = (0, -g) - (k/m)\|v\|v
	3	"""
	4	import numpy as np
	5	import lasques
	6	import matplotlib
	7	matplotlib.use("Agg")
	8	import matplotlib.pyplot as plt
	9
	10	p = lasques.params()
	11	g = 9.81
	12	v0 = float(p.get("v0", 30.0))
	13	angle_deg = float(p.get("angle_deg", 45.0))
	14	km = float(p.get("km", 0.0))
	15	dt = float(p.get("dt", 1e-3))
	16
	17	th = np.radians(angle_deg)
	18	s = np.array([0.0, 0.0, v0np.cos(th), v0np.sin(th)])
	19
	20	def deriv(s, km):
	21	x, y, vx, vy = s
	22	spd = np.hypot(vx, vy)
	23	return np.array([vx, vy, -kmspdvx, -g - kmspdvy])
	24
	25	def rk4_step(s, dt, km):
	26	k1 = deriv(s, km)
	27	k2 = deriv(s + 0.5dtk1, km)
	28	k3 = deriv(s + 0.5dtk2, km)
	29	k4 = deriv(s + dt*k3, km)
	30	return s + dt(k1 + 2k2 + 2*k3 + k4)/6.0
	31
	32	traj = [s.copy()]
	33	t = 0.0
	34	while s[1] >= 0.0:
	35	s_new = rk4_step(s, dt, km)
	36	if s_new[1] < 0.0:
	37	# linear interpolation of ground crossing
	38	y0, y1 = s[1], s_new[1]
	39	frac = y0 / (y0 - y1)
	40	s_land = s + frac*(s_new - s)
	41	s = s_land
	42	traj.append(s)
	43	t += frac*dt
	44	break
	45	s = s_new
	46	t += dt
	47	traj.append(s.copy())
	48	if t > 60:
	49	raise RuntimeError("no ground crossing")
	50
	51	traj = np.array(traj)
	52	range_ = s[0]
	53	tof = t
	54	apex = traj[:, 1].max()
	55
	56	lasques.log_metric("y", s[1], step=0)
	57	lasques.set_metrics({
	58	"range_m": range_,
	59	"time_of_flight_s": tof,
	60	"apex_m": apex,
	61	"km": km, "v0": v0, "angle_deg": angle_deg, "dt": dt,
	62	})
	63	lasques.save_json({"traj_x": traj[:,0].tolist(), "traj_y": traj[:,1].tolist()}, "trajectory")
	64	lasques.summary(f"km={km}, angle={angle_deg}, dt={dt}: range={range_:.4f} m, apex={apex:.3f} m, tof={tof:.3f} s")