diff --git a/experiment/newton/pendulum.c b/experiment/newton/pendulum.c index 066f4b21ca0905ece935c628a2b847f86db1bc5c..d8fb8ccaead0392c2aa87c9680d9d84d48aaf09d 100644 --- a/experiment/newton/pendulum.c +++ b/experiment/newton/pendulum.c @@ -35,15 +35,25 @@ point point_sub(point lhs, point rhs) return res; } +point point_neg(point lhs) +{ + return point_new(-lhs.x, -lhs.y); +} + point point_mul(double lhs, point rhs) { point res = {lhs * rhs.x, lhs * rhs.y}; return res; } +double point_scalar_product(point lhs, point rhs) +{ + return lhs.x * rhs.x + lhs.y * rhs.y; +} + double point_norm_sqr(point lhs) { - return lhs.x * lhs.x + lhs.y * lhs.y; + return point_scalar_product(lhs, lhs); } double point_norm(point lhs) @@ -51,6 +61,11 @@ double point_norm(point lhs) return sqrt(point_norm_sqr(lhs)); } +point point_normalize(point lhs) +{ + return point_mul(1.0 / point_norm(lhs), lhs); +} + point compute_deriv(double dt, point p, point p_1) { return point_mul(1.0 / dt, point_sub(p, p_1)); @@ -58,10 +73,12 @@ point compute_deriv(double dt, point p, point p_1) point compute_force(point pos, double g, double mass) { - double angle = atan2(pos.y, pos.x); + // double angle = atan2(pos.y, pos.x); double amplitude = g * mass; point grav = point_new(0.0, -amplitude); - point string = point_from_polar(amplitude, angle); + point neg = point_neg(point_normalize(pos)); + double prod = point_scalar_product(neg, grav); + point string = point_mul(-prod, neg); return point_add(grav, string); } @@ -83,8 +100,8 @@ static void render(struct gfx_context_t *context, point pos, double length) int cx = context->width / 2; int cy = context->height / 2; render_point(context, cx, cy, COLOR_RED); - int x = cx + pos.x / length * context->width / 2; - int y = cy + pos.y / length * context->height / 2; + int x = cx + pos.x / length * context->width / 4; + int y = cy - pos.y / length * context->height / 4; render_point(context, x, y, COLOR_WHITE); gfx_present(context); @@ -109,7 +126,7 @@ int main(int argc, char **argv) { mass = atof(argv[1]); length = atof(argv[2]); - debug = (0 == strcmp("-d", argv[2])); + debug = (0 == strcmp("-d", argv[3])); } else { @@ -128,11 +145,12 @@ int main(int argc, char **argv) return EXIT_FAILURE; } - double angle = (double)rand() / (double)RAND_MAX * 2.0 * M_PI; + // double angle = (double)rand() / (double)RAND_MAX * 2.0 * M_PI; + double angle = (double)-M_PI/2.0+0.1; point p = point_from_polar(length, angle); point p_1 = p; // TODO: change that.... point p_2 = p; // TODO: change that.... - double dt = 0.01; + double dt = 0.001; struct timespec time = {0, dt * 1.0e9}; @@ -154,7 +172,7 @@ int main(int argc, char **argv) point vel = compute_deriv(dt, p, p_1); point vel_1 = compute_deriv(dt, p_1, p_2); point acc = compute_deriv(dt, vel, vel_1); - if (i > 0) + if (i >= 0) { if (debug) { @@ -171,7 +189,7 @@ int main(int argc, char **argv) nanosleep(&time, NULL); point p_tmp = p; point new_acc = point_mul(1.0 / mass, compute_force(p, GRAV, mass)); - p = verlet(p, p_1, new_acc, dt); // add some noise of the order of + p = verlet(p, p_1, new_acc, dt); p_2 = p_1; p_1 = p_tmp; }