function AsteroidShowdown(sizeFactor, s)

w = s;
h = s;

%% "gravitational" constant to make for a good picture
G = .01;
dt = 1;

%% Calculate radius of a planet by index
function r = planetRadius(i)
    r = sqrt(planetMass(i));
end

%% Clear the figure, and draw our planets
function drawPlanets()
    clf('reset');
    axis([0 w 0 h]);
    %axis equal;
    hold on;
    for j = 1:numPlanets
        if planetPosition(j,1) == -1
            continue;
        end
        r = planetRadius(j);
        if r == 0
            continue;
        end
        pos = [(planetPosition(j,:)-r/2) r r];
        rectangle('position', pos, 'curvature', [1 1], ...
            'facecolor', planetColor(j,:));
    end
end

%% Rotation matrix generator
function rot = rotation(rad)
    rot = [cos(rad) -sin(rad); sin(rad) cos(rad)];
end

%% Reflect vector across a line
function vPrime = reflect(v, n)
    n = n/norm(n);
    vPrime = 2 * v + 2 * n * dot(v,n);
end

%% Move the planets due to gravity
function movePlanets()
    for j = 1:numPlanets
        if planetPosition(j,1) == -1
            continue;
        end
        Eforce = [0 0];
        for k = 1:numPlanets
            if j == k || planetPosition(k,1) == -1
                continue;
            end
            r = planetPosition(k,:) - planetPosition(j,:);
            if abs(r) < 5
                continue;
            end
            force = G * planetMass(k) / ...
                norm(r/5)^3 * (r/5);
            Eforce = Eforce + force;
        end
        planetVelocity(j,:) = planetVelocity(j,:) + Eforce;
        if j ~= 1
            planetVelocity(j,:) = [0 0];
        end
        planetPosition(j,:) = planetPosition(j,:) + planetVelocity(j,:) .* dt;
        
        if min(planetPosition(j,:)) < -50 || max(planetPosition(j,:) > s + 50)
            planetPosition(j,:) = [-1 -1];
            activePlanets = activePlanets - 1;
        end
    end
end

%% Crash planets into each other!
function checkCollisions()
    for j = 1:numPlanets
        if planetPosition(j,1) == -1
            continue;
        end
        for k = 1:numPlanets
            if planetPosition(k,1) == -1
                continue;
            end
            if j == k
                continue;
            end
            if abs(norm(planetPosition(j,:) - planetPosition(k,:))) < ...
                    (planetRadius(j) + planetRadius(k))/2
                %fprintf('Collision between %d and %d\n', j, k);
                if planetMass(j) < 1.25*planetMass(k)
                    if planetMass(k) < 1.25*planetMass(j)
                        vJ = planetVelocity(j,:);
                        vK = planetVelocity(k,:);
                        pJ = planetPosition(j,:);
                        pK = planetPosition(k,:);
                        bounceNorm = (pJ-pK) * rotation(pi/2);
                        %planetVelocity(j,:) = reflect(vJ, bounceNorm)/2;
                        %planetVelocity(k,:) = reflect(vK, bounceNorm)/2;
                        planetVelocity(j,:) = [vJ(2) -vJ(1)];
                        planetVelocity(k,:) = [vK(2) -vK(1)];
                        %fprintf('vJ from %f,%f to %f,%f\n', vJ(1), vJ(2), planetVelocity(j,1), planetVelocity(j,2));
                        continue;
                    else
                        % handled when [j,k] = [k,j]
                        continue;
                    end
                end
                Em = planetMass(j) + planetMass(k);
                pJ = planetMass(j) / Em;
                pK = planetMass(k) / Em;
                planetPosition(j,:) = pJ * planetPosition(j,:) + pK * planetPosition(k,:);
                planetPosition(k,:) = [-1 -1];
                planetMass(j) = norm([planetMass(j) planetMass(k)]);
                planetColor(j) = (planetColor(j) + planetColor(k)) / 2;
                planetVelocity(j,:) = (planetVelocity(j,:) * pJ + ...
                    planetVelocity(k,:) * pK) / 2;
                activePlanets = activePlanets - 1;
            end
        end
    end
end

%% Pick initial conditions
i = 0;
while true
    i = i + 1;
    input = ginput(1);
    if size(input,1) == 0
        break;
    end
    numPlanets = i;
    planetMass(i) = (rand() * 200 + 50) * sizeFactor;
    r = planetRadius(i);
    range = [w h] - 2*r;
    %planetPosition(i,:) = [rand() rand()] .* range + r;
    planetPosition(i,:) = input;
    planetVelocity(i,:) = [0 0]; %([planetPosition(i,2) planetPosition(i,1)] - [s/2 s/2]) ./ planetMass(i);
    planetColor(i,:) = [rand() rand() rand()];
    drawPlanets();
end
activePlanets = numPlanets
planetVelocity(1,:) = (ginput(1) - planetPosition(1,:)) / 50;

%% MAIN LOOP
while planetPosition(1,1) > -1
    %% move planets
    movePlanets()
    %planetVelocity(1,:) = [0 0];
    checkCollisions()

    %% draw planets
    drawPlanets()
    drawnow;
end
end

