dt = collision_time(N, p_next, dim, gamma, h, &col_p1, &col_p2, &min_r);
        //printf("%Lg %.20g %.20g %.20g\n", dt, min_r, p[col_p1].x, p[col_p2].x);

	if(dt > 0.0)
	{
            copy_state(N, p, p_next);
            dt += dt/1000;
            if (dt < h)
            {
                hprime = h - dt;
                euler_velocity(N, p_next, dW, f, hprime, dif, ktm);
                euler_displacement(N, p_next, hprime);
            }
	    collision_action(&p_next[col_p1], &p_next[col_p2], dim);
            for(int j=0; j<N; j++)
                boundary(dim, &p_next[j]);
	}
	copy_state(N, p_next, p);
	null_force(N, f);