From 44d868b0684ef88d830f2e99624b826fc4c3507e Mon Sep 17 00:00:00 2001
From: rodri <rgl@antares-labs.eu>
Date: Tue, 16 Jan 2024 23:02:46 +0000
Subject: add air drag.

---
 dat.h  |  6 +++++-
 main.c | 38 ++++++++++++++++++++++++++------------
 2 files changed, 31 insertions(+), 13 deletions(-)

diff --git a/dat.h b/dat.h
index 24799b8..da2038d 100644
--- a/dat.h
+++ b/dat.h
@@ -1,4 +1,6 @@
-#define Eg	9.81
+#define Eg	9.81	/* earth's gravity (m·s⁻²) */
+#define Cd	0.45	/* drag coefficient for a sphere */
+#define ρ	1.293	/* air density (kg·m⁻³) */
 #define PIX2M	0.001
 #define M2PIX	(1.0/PIX2M)
 
@@ -6,6 +8,7 @@ enum {
 	Stheta = 0,
 	Spos,
 	Svel,
+	Sdrag,
 	Sdeltax,
 	Seta,
 	SLEN,
@@ -17,4 +20,5 @@ struct Projectile
 {
 	Point2 p, v;
 	double mass;
+	double r;
 };
diff --git a/main.c b/main.c
index 1547898..42fbaf8 100644
--- a/main.c
+++ b/main.c
@@ -15,6 +15,7 @@ RFrame worldrf;
 Projectile ball;
 double t0, Δt;
 double v0;
+double A;	/* area of the projectile that meets the air */
 Point2 target;
 
 char stats[SLEN][64];
@@ -46,12 +47,28 @@ fromscreen(Point p)
 	return rframexform(Pt2(p.x, p.y, 1), worldrf);
 }
 
+void
+∫(double Δt)
+{
+	Point2 Fd;	/* drag force */
+
+	Fd = mulpt2(mulpt2(ball.v, -vec2len(ball.v)), 0.5 * Cd*ρ*A);	/* ½CdρAv² */
+	ball.v = addpt2(ball.v, mulpt2(addpt2(Vec2(0,-Eg), divpt2(Fd, ball.mass)), Δt));
+	ball.p = addpt2(ball.p, mulpt2(ball.v, Δt));
+	snprint(stats[Spos], sizeof(stats[Spos]), "p: %v", ball.p);
+	snprint(stats[Sdrag], sizeof(stats[Sdrag]), "Fd: %v", Fd);
+	if(ball.p.y <= (2+1)*M2PIX){
+		ball.p.y = (2+1)*M2PIX;
+		ball.v = Vec2(0,0);
+	}
+}
+
 void
 drawstats(void)
 {
 	int i;
 
-	snprint(stats[Svel], sizeof(stats[Svel]), "v: %gm/s", vec2len(ball.v));
+	snprint(stats[Svel], sizeof(stats[Svel]), "|v|: %gm/s", vec2len(ball.v));
 	snprint(stats[Sdeltax], sizeof(stats[Sdeltax]), "Δx: %gm", target.x-ball.p.x);
 	for(i = 0; i < nelem(stats); i++)
 		stringn(screen, addpt(screen->r.min, Pt(10, font->height*i+1)), statc, ZP, font, stats[i], sizeof(stats[i]));
@@ -86,7 +103,8 @@ mmb(void)
 	switch(menuhit(2, mc, &menu, nil)){
 	case SETV0:
 		enter("v0(m/s):", buf, sizeof(buf), mc, kc, nil);
-		v0 = strtod(buf, 0);
+		if(buf[0] != 0)
+			v0 = strtod(buf, nil);
 		break;
 	}
 }
@@ -124,7 +142,7 @@ mouse(void)
 	if(ball.p.y <= (2+1)*M2PIX){
 		p = subpt2(fromscreen(mc->xy), ball.p);
 		θ = atan2(p.y, p.x);
-		snprint(stats[Stheta], sizeof(stats[Stheta]), "θ: %g°", θ*180/PI);
+		snprint(stats[Stheta], sizeof(stats[Stheta]), "θ: %g°", θ/DEG);
 		dist = v0*v0*sin(2*θ)/Eg;
 		target = Pt2(ball.p.x+dist, 0, 1);
 		eta = 2*v0*sin(θ)/Eg;
@@ -204,6 +222,8 @@ threadmain(int argc, char *argv[])
 	ball.p = Pt2((2+1)*M2PIX,(2+1)*M2PIX,1);
 	ball.v = Vec2(0, 0);
 	ball.mass = 106000;
+	ball.r = 0.1;		/* 10cm */
+	A = 2*PI*ball.r*ball.r;	/* ½(4πr²) */
 	v0 = 1640; /* Paris Gun's specs */
 
 	scrsync = chancreate(1, 0);
@@ -227,14 +247,8 @@ threadmain(int argc, char *argv[])
 		case 2: key(r); break;
 		case 3: redraw(); break;
 		}
-		Δt = (nsec()-t0)/1e9;
-		ball.v = addpt2(ball.v, mulpt2(Vec2(0,-Eg), Δt));
-		ball.p = addpt2(ball.p, mulpt2(ball.v, Δt));
-		snprint(stats[Spos], sizeof(stats[Spos]), "p: %v", ball.p);
-		if(ball.p.y <= (2+1)*M2PIX){
-			ball.p.y = (2+1)*M2PIX;
-			ball.v = Vec2(0,0);
-		}
-		t0 += Δt*1e9;
+		Δt = (nsec()-t0);
+		∫(Δt/1e9);
+		t0 += Δt;
 	}
 }
-- 
cgit v1.2.3