if (!((rt[0] <= EPSILON) && (uveloc.dot(norm) > EPSILON))) {
                                normal.set(norm);
                                point.set(Nc);
                                lamda[0] = rt4;
                                BallNr = i;
                            }
                    }
                }

                if (TestIntersionCylinder(cyl3, oldPos[i], uveloc, rt, norm, Nc)) {
                    rt4 = rt[0] * RestTime[0] / rt2;
                    if (rt4 <= lamda[0]) {
                        if (rt4 <= RestTime[0] + EPSILON)
                            if (!((rt[0] <= EPSILON) && (uveloc.dot(norm) > EPSILON))) {
                                normal.set(norm);
                                point.set(Nc);
                                lamda[0] = rt4;
                                BallNr = i;
                            }
                    }
                }
            }

            //After all balls were tested with planes/cylinders test for collision
            //between them and replace if collision time smaller
            if (findBallCol(Pos2, BallTime, RestTime, BallColNr1, BallColNr2) == 1) {
                if (soundsEnabled)
                    audioSample.play();

                if ((lamda[0] == 10000) || (lamda[0] > BallTime[0])) {
                    RestTime[0] = RestTime[0] - BallTime[0];
                    Tuple3d pb1 = new Tuple3d(),
                            pb2 = new Tuple3d(),
                            xaxis = new Tuple3d(),
                            U1x = new Tuple3d(),
                            U1y = new Tuple3d(),
                            U2x = new Tuple3d(),
                            U2y = new Tuple3d(),
                            V1x = new Tuple3d(),
                            V1y = new Tuple3d(),
                            V2x = new Tuple3d(),
                            V2y = new Tuple3d();
                    double a,b;

                    pb1.scaleAdd(BallTime[0], arrayVel[BallColNr1[0]], oldPos[BallColNr1[0]]);
                    pb2.scaleAdd(BallTime[0], arrayVel[BallColNr2[0]], oldPos[BallColNr2[0]]);
                    xaxis.sub(pb2, pb1);
                    xaxis.normalize();

                    a = xaxis.dot(arrayVel[BallColNr1[0]]);
                    U1x.scaleAdd(a, xaxis);
                    U1y.sub(arrayVel[BallColNr1[0]], U1x);

                    xaxis.sub(pb1, pb2);
                    xaxis.normalize();

                    b = xaxis.dot(arrayVel[BallColNr2[0]]);
                    U2x.scaleAdd(b, xaxis);
                    U2y.sub(arrayVel[BallColNr2[0]], U2x);

                    V1x.add(U1x, U2x);
                    V1x.sub(new Tuple3d(U1x.x - U2x.x, U1x.y - U2x.y, U1x.z - U2x.z));
                    V1x.scale(.5);

                    V2x.add(U1x, U2x);
                    V2x.sub(new Tuple3d(U2x.x - U1x.x, U2x.y - U1x.y, U2x.z - U1x.z));
                    V2x.scale(.5);

                    V1y.set(U1y);
                    V2y.set(U2y);

                    for (int j = 0; j < nrOfBalls; j++)
                        arrayPos[j].scaleAdd(BallTime[0], arrayVel[j], oldPos[j]);

                    arrayVel[BallColNr1[0]].add(V1x, V1y);
                    arrayVel[BallColNr2[0]].add(V2x, V2y);

                    //Update explosion array
                    for (int j = 0; j < 20; j++) {
                        if (explosions[j].alpha <= 0) {
                            explosions[j].alpha = 1;
                            explosions[j].position = arrayPos[BallColNr1[0]];
                            explosions[j].scale = 1;
                            break;
                        }
                    }
                    continue;
                }
            }
            //End of tests
            //If test occured move simulation for the correct timestep
            //and compute response for the colliding ball
            if (lamda[0] != 10000) {
                RestTime[0] -= lamda[0];

                for (int j = 0; j < nrOfBalls; j++)
                    arrayPos[j].scaleAdd(lamda[0], arrayVel[j], oldPos[j]);

                rt2 = arrayVel[BallNr].length();
                arrayVel[BallNr].normalize();

                normal.scale(-2 * normal.dot(arrayVel[BallNr]));
                arrayVel[BallNr].add(normal, arrayVel[BallNr]);
                arrayVel[BallNr].normalize();
                arrayVel[BallNr].scale(rt2);

                // Update explosion array
                for (int j = 0; j < 20; j++) {
                    if (explosions[j].alpha <= 0) {
                        explosions[j].alpha = 1;
                        explosions[j].position = point;
                        explosions[j].scale = 1;
                        break;
                    }
                }
            } else
                RestTime[0] = 0;
        }
    }

    public void reshape(GLAutoDrawable drawable,
                        int xstart,
                        int ystart,
                        int width,
                        int height) {
        GL gl = drawable.getGL();

        height = (height == 0) ? 1 : height;

        gl.glViewport(0, 0, width, height);         // Reset The Current Viewport
        gl.glMatrixMode(GL.GL_PROJECTION);       // Select The Projection Matrix
        gl.glLoadIdentity();                     // Reset The Projection Matrix

        // Calculate The Aspect Ratio Of The Window
        glu.gluPerspective(50.0f, (float) width / height, 10.f, 1700.0f);

        gl.glMatrixMode(GL.GL_MODELVIEW);        // Select The Modelview Matrix
        gl.glLoadIdentity();
    }

    public void displayChanged(GLAutoDrawable drawable,
                               boolean modeChanged,
                               boolean deviceChanged) {
    }

    /*
     * Intersetion tests
     */

    /**
     * Fast Intersection Function between ray/plane
     */
    private boolean TestIntersionPlane(Plane plane, Tuple3d position, Tuple3d direction,
                                       double[] lamda, Tuple3d pNormal) {

        double dotProduct = direction.dot(plane.normal);
        double l2;

        //determine if ray paralle to plane
        if ((dotProduct < EPSILON) && (dotProduct > -EPSILON))
            return false;

        Tuple3d substract = new Tuple3d(plane.position);
        substract.sub(position);
        l2 = (plane.normal.dot(substract)) / dotProduct;

        if (l2 < -EPSILON)
            return false;

        pNormal.set(plane.normal);
        lamda[0] = l2;
        return true;
    }

    /**
     * Fast Intersection Function between ray/cylinder
     */
    private boolean TestIntersionCylinder(Cylinder cylinder, Tuple3d position,
                                          Tuple3d direction, double[] lamda,
                                          Tuple3d pNormal, Tuple3d newposition) {
        Tuple3d RC = new Tuple3d(),
                HB = new Tuple3d(),
                n = new Tuple3d(),
                O = new Tuple3d();
        double d,t,s,
                ln,in,out;

        RC.sub(position, cylinder.position);
        n.cross(direction, cylinder.axis);

        ln = n.length();

        if ((ln < EPSILON) && (ln > -EPSILON))
            return false;

        n.normalize();
        d = Math.abs(RC.dot(n));

        if (d <= cylinder.radius) {

            O.cross(RC, cylinder.axis);
            t = -O.dot(n) / ln;
            O.cross(n, cylinder.axis);
            O.normalize();
            s = Math.abs(Math.sqrt(cylinder.radius * cylinder.radius - d * d) / direction.dot(O));

            in = t - s;
            out = t + s;

            if (in < -EPSILON) {
                if (out < -EPSILON)
                    return false;
                else
                    lamda[0] = out;
            } else if (out < -EPSILON) {
                lamda[0] = in;
            } else if (in < out)
                lamda[0] = in;
            else
                lamda[0] = out;

            newposition.scaleAdd(lamda[0], direction, position);
            HB.sub(newposition, cylinder.position);
            pNormal.scaleAdd(-HB.dot(cylinder.axis), cylinder.axis, HB);
            pNormal.normalize();
            return true;
        }
        return false;
    }

    /**
     * Find if any of the current balls intersect with each other in the current timestep.
     * @return the index of the 2 itersecting balls, the point and time of intersection
     */
    private int findBallCol(Tuple3d point, double[] TimePoint,
                            double[] Time2, int[] BallNr1, int[] BallNr2) {

        Tuple3d RelativeVClone = new Tuple3d(),
                RelativeV = new Tuple3d(),
                posi = new Tuple3d();

        double Timedummy = 10000,
                MyTime = 0,
                Add = Time2[0] / 150;
        TRay rays;

        //Test all balls against eachother in 150 small steps
        for (int i = 0; i < nrOfBalls - 1; i++) {
            for (int j = i + 1; j < nrOfBalls; j++) {

                RelativeV.sub(arrayVel[i], arrayVel[j]);
                RelativeVClone.set(RelativeV);
                RelativeVClone.normalize();
                rays = new TRay(oldPos[i], RelativeVClone);
                MyTime = 0;

                if ((rays.dist(oldPos[j])) > 40)
                    continue;

                while (MyTime < Time2[0]) {

                    MyTime += Add;
                    posi.scaleAdd(MyTime, RelativeV, oldPos[i]);
                    if (posi.distance(oldPos[j]) <= 40) {
                        point.set(posi);
                        if (Timedummy > (MyTime - Add))
                            Timedummy = MyTime - Add;

                        BallNr1[0] = i;
                        BallNr2[0] = j;
                        break;
                    }
                }
            }
        }

        if (Timedummy != 10000) {
            TimePoint[0] = Timedummy;
            return 1;
        }
        return 0;
    }

    private static class Explosion {
        public Tuple3d position = new Tuple3d();
        public float alpha;
        public float scale;
    }

    private static class Cylinder {
        public Tuple3d position;
        public Tuple3d axis;
        public double radius;
    }

    private static class Plane {
        public Tuple3d position;
        public Tuple3d normal;
    }

    private static class TRay {
        public Tuple3d p;                     // Any point on the line
        public Tuple3d v;                     // Direction of the line

        public TRay(Tuple3d p,
                    Tuple3d v) {
            this.p = new Tuple3d(p);
            this.v = new Tuple3d(v);
        }

        public double dist(Tuple3d point) {
            double lambda = v.dot(new Tuple3d(point.x - p.x,
                    point.y - p.y,
                    point.z - p.z));
            Tuple3d point2 = new Tuple3d();
            point2.scaleAdd(lambda, v, p);
            return point.distance(point2);
        }
    }
}