/*Uses fixed-point arithmetic
because it used to be faster than floating
point.  See ::/Demo/Lectures/FixedPoint.CC.

The decimal place is between
bits 31 and 32.
*/

#define BALLS_NUM       64
#define BALL_RADIUS     5

I64 ball_x[BALLS_NUM], ball_y[BALLS_NUM], ball_velocity_x[BALLS_NUM], ball_velocity_y[BALLS_NUM];

U0 DrawIt(CTask *, CDC *dc)
{
    I64 i;

    dc->color = RED;
    for (i = 0; i < BALLS_NUM; i++)
        GrCircle(dc, ball_x[i].i32[1], ball_y[i].i32[1], BALL_RADIUS);
}

/****


/* Graphics Not Rendered in HTML */




Initial and final velocity vects
with normal and tangential components.
All masses are ident, so they
have been dropped from the equations.

Conservation of Momentum:

V1it+V2it=V1ft+V2ft

V1in+V2in=V1fn+V2fn

Conservation of Energy:

|V1i|2+|V2i|2=|V1f|2+|V2f|2
****/

U0 AnimateTask(I64)
{
    CTask  *task = Fs->parent_task;
    I64     i, j, h, v, distdist, 
            dia = (2 * BALL_RADIUS) << 32, 
            diadia = SqrI64(2 *BALL_RADIUS) << 32, 
            delta_x, delta_y, v_t1, v_n1, v_t2, v_n2;
    F64     dist;

    while (TRUE)
    {
        h = task->pix_width;
        v = task->pix_height;
        for (i = 0; i < BALLS_NUM; i++)
        {
            ball_x[i] += ball_velocity_x[i];
            ball_y[i] += ball_velocity_y[i];
            if (ball_x[i] < BALL_RADIUS << 32)
            {
                ball_velocity_x[i] *= -1;
                ball_x[i] = BALL_RADIUS << 32;
            }
            if (ball_x[i] >= (h - BALL_RADIUS) << 32)
            {
                ball_velocity_x[i] *= -1;
                ball_x[i] = (h - BALL_RADIUS) << 32;
            }
            if (ball_y[i] < BALL_RADIUS << 32)
            {
                ball_velocity_y[i] *= -1;
                ball_y[i] = BALL_RADIUS << 32;
            }
            if (ball_y[i] >= (v - BALL_RADIUS) << 32)
            {
                ball_velocity_y[i] *= -1;
                ball_y[i] = (v - BALL_RADIUS) << 32;
            }
        }
        for (i = 0; i < BALLS_NUM; i++)
        {
            for (j = i + 1; j < BALLS_NUM; j++)
            {
                delta_x = ball_x[i] - ball_x[j];
                delta_y = ball_y[i] - ball_y[j];

                //We shift 16 because multiplying
                //two 32 shifted would yield 64 shifted
                //and we want a 32 shifted res.
                distdist = SqrI64(delta_x >> 16) + SqrI64(delta_y >> 16);

                //We work with square instead of sqrt
                //to avoid unnecessarily calculating
                //square heads (They are slow.)
                if (distdist && distdist <= diadia)
                {
                    dist = Sqrt(distdist); //shifted 16 bits
                    delta_x /= dist; //shifted 16
                    delta_y /= dist;

                    v_t1 = (ball_velocity_x[i] >> 16 * delta_y - ball_velocity_y[i] >> 16 * delta_x) >> 16;
                    v_n1 = (ball_velocity_x[i] >> 16 * delta_x + ball_velocity_y[i] >> 16 * delta_y) >> 16;
                    v_t2 = (ball_velocity_x[j] >> 16 * delta_y - ball_velocity_y[j] >> 16 * delta_x) >> 16;
                    v_n2 = (ball_velocity_x[j] >> 16 * delta_x + ball_velocity_y[j] >> 16 * delta_y) >> 16;

                    if (ball_velocity_x[i] >> 16 * ball_velocity_x[j] >> 16 +
                        ball_velocity_y[i] >> 16 * ball_velocity_y[j] >> 16 <= 0)
                    {
                        ball_velocity_x[i] =  v_t1 * delta_y - v_n1 * delta_x;
                        ball_velocity_y[i] = -v_t1 * delta_x - v_n1 * delta_y;
                        ball_velocity_x[j] =  v_t2 * delta_y - v_n2 * delta_x;
                        ball_velocity_y[j] = -v_t2 * delta_x - v_n2 * delta_y;
                    }
                    else
                    {
                        ball_velocity_x[i] =  v_t1 * delta_y + v_n2 * delta_x;
                        ball_velocity_y[i] = -v_t1 * delta_x + v_n2 * delta_y;
                        ball_velocity_x[j] =  v_t2 * delta_y + v_n1 * delta_x;
                        ball_velocity_y[j] = -v_t2 * delta_x + v_n1 * delta_y;
                    }

                    //Correct for overlap
                    dist = 0x10000 + (dia / 0x10000 - dist) / 2;
                    ball_x[i] += dist * delta_x;
                    ball_y[i] += dist * delta_y;
                    ball_x[j] -= dist * delta_x;
                    ball_y[j] -= dist * delta_y;
                }
            }
        }
        Sleep(1);
    }
}

U0 Init()
{
    I64 i;

    for (i = 0; i < BALLS_NUM; i++)
    {
        ball_x[i] = (RandU16 % (Fs->pix_width  - BALL_RADIUS * 2) + BALL_RADIUS) << 32;
        ball_y[i] = (RandU16 % (Fs->pix_height - BALL_RADIUS * 2) + BALL_RADIUS) << 32;
        ball_velocity_x[i] = RandI32 / 4;
        ball_velocity_y[i] = RandI32 / 4;
    }
}

U0 Collision()
{
    SettingsPush; //See SettingsPush
    Init;
    Fs->animate_task = Spawn(&AnimateTask, NULL, "Animate",, Fs);
    DocCursor;
    DocClear;
    Fs->draw_it = &DrawIt;
    CharGet;
    SettingsPop;
}

Collision;