diff --git a/Makefile b/Makefile
index 8bde1d4..541b388 100644
--- a/Makefile
+++ b/Makefile
@@ -13,13 +13,13 @@ LIBS=$(GTK_LIBS) -lm
 PROGS=balls
 OBJS=balls.o c_index.o game.o gravity.o spaceship.o main.o
 
-# dependencies (gcc -MM *.c)
-balls.o: balls.cc game.h balls.h gravity.h
-c_index.o: c_index.cc balls.h game.h c_index.h
+# dependencies (gcc -MM *.cc)
+balls.o: balls.cc game.h balls.h vec2d.h gravity.h
+c_index.o: c_index.cc balls.h vec2d.h game.h c_index.h
 game.o: game.cc game.h
-gravity.o: gravity.cc gravity.h game.h
-main.o: main.cc game.h balls.h c_index.h gravity.h spaceship.h
-spaceship.o: spaceship.cc balls.h game.h
+gravity.o: gravity.cc gravity.h balls.h vec2d.h game.h
+main.o: main.cc game.h balls.h vec2d.h c_index.h gravity.h spaceship.h
+spaceship.o: spaceship.cc balls.h vec2d.h game.h
 stats.o: stats.cc
 
 .PHONY: run
diff --git a/balls.cc b/balls.cc
index a514d1b..8892e50 100644
--- a/balls.cc
+++ b/balls.cc
@@ -18,13 +18,19 @@ unsigned int v_angle_max = 100;
 ball * balls = nullptr;
 unsigned int n_balls = 50;
 
-static void random_velocity(ball * p) {
+static vec2d random_velocity() {
     double r2;
+    vec2d v;
     do {
-	p->v_x = v_min + rand() % (v_max + 1 - v_min);
-	p->v_y = v_min + rand() % (v_max + 1 - v_min);
-	r2 = p->v_x*p->v_x + p->v_y*p->v_y;
+	v.x = v_min + rand() % (v_max + 1 - v_min);
+	v.y = v_min + rand() % (v_max + 1 - v_min);
+	r2 = vec2d::dot(v,v);
     } while (r2 > v_max*v_max || r2 < v_min*v_min);
+    if (rand() % 2)
+	v.x = -v.x;
+    if (rand() % 2)
+	    v.y = -v.y;
+    return v;
 }
 
 void balls_init_state () {
@@ -33,13 +39,9 @@ void balls_init_state () {
     int h = height < 2*border ? 1 : height - 2*border;
 
     for (unsigned int i = 0; i < n_balls; ++i) {
-	balls[i].x = border + rand() % w;
-	balls[i].y = border + rand() % h;
-	random_velocity(balls + i);
-	if (rand() % 2)
-	    balls[i].v_x = -balls[i].v_x;
-	if (rand() % 2)
-	    balls[i].v_y = -balls[i].v_y;
+	balls[i].position.x = border + rand() % w;
+	balls[i].position.y = border + rand() % h;
+	balls[i].velocity = random_velocity();
 	balls[i].radius = radius_min + rand() % (radius_max + 1 - radius_min);
 	unsigned int v_angle_360 = (v_angle_min + rand() % (v_angle_max + 1 - v_angle_min)) % 360;
 	balls[i].v_angle = 2*M_PI*v_angle_360/360;
@@ -48,36 +50,32 @@ void balls_init_state () {
 }
 
 void ball_update_state (ball * p) {
-    struct gravity_vector g;
-    gravity_get_vector (&g, p);
+    vec2d g = gravity_vector (p);
 
-    p->x += delta*p->v_x + delta*delta*g.x/2.0;
-    p->v_x += delta*g.x;
+    p->position += delta*p->velocity + delta*delta*g/2.0;
+    p->velocity += delta*g;
 
-    p->y += delta*p->v_y + delta*delta*g.y/2.0;
-    p->v_y += delta*g.y;
-
-    if (p->x + p->radius > width) { /* right wall */
-	if (p->v_x > 0) {
-	    p->x -= p->x + p->radius - width;
-	    p->v_x = -p->v_x;
+    if (p->position.x + p->radius > width) { /* right wall */
+	if (p->velocity.x > 0) {
+	    p->position.x -= p->position.x + p->radius - width;
+	    p->velocity.x = -p->velocity.x;
 	}
-    } else if (p->x < p->radius) { /* left wall */
-	if (p->v_x < 0) {
-	    p->x += p->radius - p->x;
-	    p->v_x = -p->v_x;
+    } else if (p->position.x < p->radius) { /* left wall */
+	if (p->velocity.x < 0) {
+	    p->position.x += p->radius - p->position.x;
+	    p->velocity.x = -p->velocity.x;
 	}
     } 
 
-    if (p->y + p->radius > height) { /* bottom wall */
-	if (p->v_y > 0) {
-	    p->y -= p->y + p->radius - height;
-	    p->v_y = -p->v_y;
+    if (p->position.y + p->radius > height) { /* bottom wall */
+	if (p->velocity.y > 0) {
+	    p->position.y -= p->position.y + p->radius - height;
+	    p->velocity.y = -p->velocity.y;
 	}
-    } else if (p->y < p->radius) { /* top wall */
-	if (p->v_y < 0) {
-	    p->y += p->radius - p->y;
-	    p->v_y = -p->v_y;
+    } else if (p->position.y < p->radius) { /* top wall */
+	if (p->velocity.y < 0) {
+	    p->position.y += p->radius - p->position.y;
+	    p->velocity.y = -p->velocity.y;
 	}
     } 
     p->angle += delta*p->v_angle;
@@ -88,39 +86,34 @@ void ball_update_state (ball * p) {
 }
 
 void ball_elastic_collision (ball * p, ball * q) {
-    double dx = q->x - p->x;
-    double dy = q->y - p->y;
-    double d2 = dx*dx + dy*dy;
+    vec2d dp = q->position - p->position;
+    double d2 = vec2d::dot(dp,dp);
     double r = p->radius + q->radius;
     if (d2 <= r*r) {
-	double dv_x = q->v_x - p->v_x;
-	double dv_y = q->v_y - p->v_y;
+	vec2d dv = q->velocity - p->velocity;
 
 	double mp = p->radius * p->radius;
 	double mq = q->radius * q->radius;
 
-	double f = dv_x*dx + dv_y*dy;
+	double f = vec2d::dot(dv,dp);
 
 	if (f < 0) {
 	    f /= d2*(mp + mq);
-	    p->v_x += 2*mq*f*dx;
-	    p->v_y += 2*mq*f*dy;
-
-	    q->v_x -= 2*mp*f*dx;
-	    q->v_y -= 2*mp*f*dy;
+	    p->velocity += 2*mq*f*dp;
+	    q->velocity -= 2*mp*f*dp;
 	}
     }
 }
 
 void ball_reposition (ball * b) {
-    if (b->x < b->radius)
-	b->x = b->radius;
-    else if (b->x + b->radius > width)
-	b->x = width - b->radius;
-    if (b->y < b->radius)
-	b->y = b->radius;
-    else if (b->y + b->radius > height)
-	b->y = height - b->radius;
+    if (b->position.x < b->radius)
+	b->position.x = b->radius;
+    else if (b->position.x + b->radius > width)
+	b->position.x = width - b->radius;
+    if (b->position.y < b->radius)
+	b->position.y = b->radius;
+    else if (b->position.y + b->radius > height)
+	b->position.y = height - b->radius;
 }
 
 const char * face_filename = 0;
@@ -234,7 +227,7 @@ static void balls_init_faces () {
 
 void ball::draw (cairo_t * cr) const {
     cairo_save (cr);
-    cairo_translate (cr, x - radius, y - radius);
+    cairo_translate (cr, position.x - radius, position.y - radius);
     cairo_set_source_surface(cr, face->get_surface (angle), 0, 0);
     cairo_paint(cr);
     cairo_restore(cr);
diff --git a/balls.h b/balls.h
index a8015b0..85a4939 100644
--- a/balls.h
+++ b/balls.h
@@ -3,16 +3,15 @@
 
 #include <gtk/gtk.h>
 
+#include "vec2d.h"
+
 class ball_face;
 
 class ball {
 public:
-    double x;
-    double y;
     unsigned int radius;
-
-    double v_x;
-    double v_y;
+    vec2d position;
+    vec2d velocity;
 
     double angle;
     double v_angle;
diff --git a/c_index.cc b/c_index.cc
index 44678e5..7f03852 100644
--- a/c_index.cc
+++ b/c_index.cc
@@ -24,24 +24,24 @@ struct bt_node * c_index = 0;
 
 static struct bt_node * c_index_init_node(struct bt_node * n, struct ball * b) {
     n->ball = b;
-    n->r.min_x = b->x - b->radius;
-    n->r.min_y = b->y - b->radius;
-    n->r.max_x = b->x + b->radius;
-    n->r.max_y = b->y + b->radius;
+    n->r.min_x = b->velocity.x - b->radius;
+    n->r.min_y = b->velocity.y - b->radius;
+    n->r.max_x = b->velocity.x + b->radius;
+    n->r.max_y = b->velocity.y + b->radius;
     n->left = 0;
     n->right = 0;
     return n;
 }
 
 static void c_index_add_ball(struct bt_node * n, const struct ball * b) {
-    if (n->r.min_x > b->x - b->radius)
-	n->r.min_x = b->x - b->radius;
-    if (n->r.min_y > b->y - b->radius)
-	n->r.min_y = b->y - b->radius;
-    if (n->r.max_x < b->x + b->radius)
-	n->r.max_x = b->x + b->radius;
-    if (n->r.max_y < b->y + b->radius)
-	n->r.max_y = b->y + b->radius;
+    if (n->r.min_x > b->velocity.x - b->radius)
+	n->r.min_x = b->velocity.x - b->radius;
+    if (n->r.min_y > b->velocity.y - b->radius)
+	n->r.min_y = b->velocity.y - b->radius;
+    if (n->r.max_x < b->velocity.x + b->radius)
+	n->r.max_x = b->velocity.x + b->radius;
+    if (n->r.max_y < b->velocity.y + b->radius)
+	n->r.max_y = b->velocity.y + b->radius;
 }
 
 static void c_index_insert(struct bt_node * t, struct bt_node * n, struct ball * b) {
@@ -53,9 +53,9 @@ static void c_index_insert(struct bt_node * t, struct bt_node * n, struct ball *
     for (;;) {
 	c_index_add_ball(t, b);
 	if (w > h) { /* horizontal split */
-	    if (b->x <= t->ball->x) {
+	    if (b->velocity.x <= t->ball->velocity.x) {
 		if (t->left) {
-		    w = t->ball->x - ref_x;
+		    w = t->ball->velocity.x - ref_x;
 		    t = t->left;
 		} else {
 		    t->left = n;
@@ -63,8 +63,8 @@ static void c_index_insert(struct bt_node * t, struct bt_node * n, struct ball *
 		}
 	    } else {
 		if (t->right) {
-		    w -= t->ball->x - ref_x;
-		    ref_x = t->ball->x;
+		    w -= t->ball->velocity.x - ref_x;
+		    ref_x = t->ball->velocity.x;
 		    t = t->right;
 		} else {
 		    t->right = n;
@@ -72,9 +72,9 @@ static void c_index_insert(struct bt_node * t, struct bt_node * n, struct ball *
 		}
 	    }
 	} else {		/* vertical split */
-	    if (b->y <= t->ball->y) {
+	    if (b->velocity.y <= t->ball->velocity.y) {
 		if (t->left) {
-		    h = t->ball->y - ref_y;
+		    h = t->ball->velocity.y - ref_y;
 		    t = t->left;
 		} else {
 		    t->left = n;
@@ -82,8 +82,8 @@ static void c_index_insert(struct bt_node * t, struct bt_node * n, struct ball *
 		}
 	    } else {
 		if (t->right) {
-		    h -= t->ball->y - ref_y;
-		    ref_y = t->ball->y;
+		    h -= t->ball->velocity.y - ref_y;
+		    ref_y = t->ball->velocity.y;
 		    t = t->right;
 		} else {
 		    t->right = n;
@@ -119,10 +119,10 @@ static struct bt_node * c_index_stack_pop() {
 }
 
 static int c_index_ball_in_rectangle(const struct bt_node * n, const struct ball * b) {
-    return n->r.min_x <= b->x + b->radius
-	&& n->r.max_x >= b->x - b->radius
-	&& n->r.min_y <= b->y + b->radius
-	&& n->r.max_y >= b->y - b->radius;
+    return n->r.min_x <= b->velocity.x + b->radius
+	&& n->r.max_x >= b->velocity.x - b->radius
+	&& n->r.min_y <= b->velocity.y + b->radius
+	&& n->r.max_y >= b->velocity.y - b->radius;
 }
 
 static int c_index_must_check(const struct bt_node * n, const struct ball * b) {
diff --git a/gravity.cc b/gravity.cc
index 9658c97..719d72b 100644
--- a/gravity.cc
+++ b/gravity.cc
@@ -4,8 +4,7 @@
 #include "gravity.h"
 #include "game.h"
 
-static double g_y = 20;
-static double g_x = 0;
+static vec2d g{.x = 0, .y = 20};
 
 static double g_r = 50;
 static double g_g = 10000000;
@@ -17,8 +16,8 @@ static int gravity_vector_init = 300;
 
 void gravity_constant_field (double x, double y) {
     constant_field = 1;
-    g_x = x;
-    g_y = y;
+    g.x = x;
+    g.y = y;
 }
 
 void gravity_newton_field (double r, double g) {
@@ -33,11 +32,11 @@ void gravity_draw (cairo_t * cr) {
 	    cairo_save(cr);
 	    cairo_new_path(cr);
 	    cairo_move_to(cr, width/2, height/2);
-	    cairo_line_to(cr, width/2 + g_x, height/2 + g_y);
+	    cairo_line_to(cr, width/2 + g.x, height/2 + g.y);
 	    cairo_set_source_rgb(cr, 1.0, 1.0, 1.0);
 	    cairo_set_line_width(cr, 1.0);
 	    cairo_stroke(cr);
-	    cairo_arc(cr, width/2 + g_x, height/2 + g_y, 3, 0, 2*M_PI);
+	    cairo_arc(cr, width/2 + g.x, height/2 + g.y, 3, 0, 2*M_PI);
 	    cairo_fill(cr);
 	    if (gravity_vector_countdown > 0)
 		--gravity_vector_countdown;
@@ -61,8 +60,8 @@ void gravity_show () {
 
 void gravity_change (double dx, double dy) {
     if (constant_field) {
-	g_x += dx;
-	g_y += dy;
+	g.x += dx;
+	g.y += dy;
 	gravity_show ();
     } else {
 	g_r += dx;
@@ -70,43 +69,32 @@ void gravity_change (double dx, double dy) {
     }
 }
 
-void gravity_get_vector (struct gravity_vector * v, const struct ball * b) {
+vec2d gravity_vector (const ball * b) {
     if (constant_field) {
-	v->x = g_x;
-	v->y = g_y;
+	return g;
     } else {
-	double dx = width/2 - b->x;
-	double dy = height/2 - b->y;
-	double r2 = dx*dx+dy*dy;
+	vec2d b_c = vec2d{width/2.0,height/2.0} - b->position;
+	double r2 = vec2d::dot(b_c,b_c);
 	if (r2 < g_r*g_r) {
-	    v->x = 0;
-	    v->y = 0;
+	    return vec2d{0,0};
 	} else {
-	    double r = sqrt(r2);
-	    v->x = g_g/r2/r*dx;
-	    v->y = g_g/r2/r*dy;
+	    return g_g/r2/sqrt(r2)*b_c;
 	}
     }
 }
 
-void gravity_collisions (struct ball * begin, struct ball * end) {
+void gravity_collisions (ball * begin, ball * end) {
     if (constant_field)
 	return;
-    for (struct ball * b = begin; b != end; ++b) {
-	double dx = b->x - width/2;
-	double dy = b->y - height/2;
-	double d2 = dx*dx + dy*dy;
+    for (ball * b = begin; b != end; ++b) {
+	vec2d b_c = b->position - vec2d{width/2.0,height/2.0};
+	double d2 = vec2d::dot(b_c, b_c);
 	double r = b->radius + g_r;
 	if (d2 <= r*r) {
-	    double dv_x = b->v_x;
-	    double dv_y = b->v_y;
-
-	    double f = dv_x*dx + dv_y*dy;
-
+	    double f = vec2d::dot(b->velocity, b_c);
 	    if (f < 0) {
 		f /= d2;
-		b->v_x -= 2*f*dx;
-		b->v_y -= 2*f*dy;
+		b->velocity -= 2*f*b_c;
 	    }
 	}
     }
diff --git a/gravity.h b/gravity.h
index 17a7964..226e536 100644
--- a/gravity.h
+++ b/gravity.h
@@ -3,15 +3,10 @@
 
 #include "balls.h"
 
-struct gravity_vector {
-    double x;
-    double y;
-};
-
 extern void gravity_constant_field (double x, double y);
 extern void gravity_newton_field (double r, double g);
 
-extern void gravity_get_vector (gravity_vector * v, const ball * b);
+extern vec2d gravity_vector (const ball * b);
 
 extern void gravity_draw (cairo_t * cr);
 extern void gravity_change (double dx, double dy);
diff --git a/main.cc b/main.cc
index e9108f5..5e24ee5 100644
--- a/main.cc
+++ b/main.cc
@@ -78,7 +78,7 @@ gint configure_event (GtkWidget *widget, GdkEventConfigure * event) {
     width = gtk_widget_get_allocated_width (widget);
     height = gtk_widget_get_allocated_height (widget);
 
-    for (struct ball * b = balls; b != balls + n_balls; ++b)
+    for (ball * b = balls; b != balls + n_balls; ++b)
 	ball_reposition (b);
     ball_reposition (&spaceship);
     return TRUE;
diff --git a/spaceship.cc b/spaceship.cc
index d11aa91..ca07d1d 100644
--- a/spaceship.cc
+++ b/spaceship.cc
@@ -10,24 +10,22 @@ int spaceship_thrust_countdown = 0;
 int spaceship_thrust_init = 50;
 
 void spaceship_init_state () {
-    spaceship.x = width/2;
-    spaceship.y = height/2;
+    spaceship.position.x = width/2;
+    spaceship.position.y = height/2;
     spaceship.radius = 30;
-    spaceship.v_x = 0;
-    spaceship.v_y = 0;
+    spaceship.velocity.x = 0;
+    spaceship.velocity.y = 0;
     spaceship.angle = 0;
     spaceship.v_angle = 0;
 }
 
 void spaceship_update_state () {
     if (spaceship_thrust > 0) {
-	double fx = cos(spaceship.angle)*spaceship_thrust*4.0;
-	double fy = sin(spaceship.angle)*spaceship_thrust*4.0;
+	vec2d f { cos(spaceship.angle)*spaceship_thrust*4.0,
+	          sin(spaceship.angle)*spaceship_thrust*4.0 };
 
-	spaceship.x += delta*delta*fx/2.0;
-	spaceship.v_x += delta*fx;
-	spaceship.y += delta*delta*fy/2.0;
-	spaceship.v_y += delta*fy;
+	spaceship.position += delta*delta*f/2.0;
+	spaceship.velocity += delta*f;
 	if (spaceship_thrust_countdown > 0)
 	    --spaceship_thrust_countdown;
 	else
@@ -40,7 +38,7 @@ void spaceship_draw (cairo_t * cr) {
     static const double one_over_sqrt_2 = 0.70710678118654752440;
     cairo_save(cr);
     cairo_set_source_rgba(cr, 0.0, 0.0, 1.0, 1.0);
-    cairo_translate(cr, spaceship.x, spaceship.y);
+    cairo_translate(cr, spaceship.position.x, spaceship.position.y);
     cairo_rotate(cr, spaceship.angle);
     cairo_arc(cr, 0, 0, spaceship.radius, 0, 2*M_PI);
     cairo_stroke(cr);
diff --git a/vec2d.h b/vec2d.h
new file mode 100644
index 0000000..1c6a3e8
--- /dev/null
+++ b/vec2d.h
@@ -0,0 +1,74 @@
+#ifndef VEC2D_H_INCLUDED
+#define VEC2D_H_INCLUDED
+
+#include <cmath>
+
+class vec2d {
+public:
+    double x;
+    double y;
+
+    vec2d & operator = (const vec2d & other) {
+	x = other.x;
+	y = other.y;
+	return *this;
+    }
+
+    vec2d & operator += (const vec2d & other) {
+	x += other.x;
+	y += other.y;
+	return *this;
+    }
+
+    vec2d & operator -= (const vec2d & other) {
+	x -= other.x;
+	y -= other.y;
+	return *this;
+    }
+
+    vec2d & operator *= (double l) {
+	x *= l;
+	y *= l;
+	return *this;
+    }
+
+    vec2d operator + (const vec2d & other) const {
+	return vec2d{x + other.x, y + other.y};
+    }
+
+    vec2d operator - (const vec2d & other) const {
+	return vec2d{x - other.x, y - other.y};
+    }
+
+    vec2d operator * (double l) const {
+	return vec2d{x*l, y*l};
+    }
+
+    vec2d operator / (double a) const {
+	return vec2d{x/a, y/a};
+    }
+
+    vec2d & rotate (double angle) {
+	double sin_a = sin(angle);
+	double cos_a = cos(angle);
+	double x1 = x;
+	double y1 = y;
+	x = cos_a*x1 - sin_a*y1;
+	y = sin_a*x1 + cos_a*y1;
+	return *this;
+    }
+
+    static double dot (const vec2d & a, const vec2d & b) {
+	return a.x*b.x + a.y*b.y;
+    }
+
+    static double cross (const vec2d & a, const vec2d & b) {
+	return a.x*b.x + a.y*b.y;
+    }
+};
+
+static inline vec2d operator * (double l, const vec2d & v) {
+    return vec2d{v.x*l, v.y*l};
+}
+
+#endif