% Galton Board Simulation - Balls dropping into bins to show normal distribution
% Each ball hits a series of pegs, going left or right with equal probability,
% then lands in a bin. The resulting distribution approximates a normal curve.

function galton_board()
    %% Parameters
    nBalls    = 300;   % total number of balls to drop
    nLevels   = 20;    % number of peg levels (determines number of bins = nLevels+1)
    nBins     = nLevels + 1;
    pegSpacing = 1;    % horizontal spacing between pegs
    levelHeight = 1;   % vertical spacing between peg levels

    % Derived geometry
    boardTop    = (nLevels + 2) * levelHeight;
    boardBottom = 0;
    binWidth    = pegSpacing;
    boardLeft   = -nBins/2 * binWidth;
    boardRight  =  nBins/2 * binWidth;

    %% Set up figure
    fig = figure('Name','Galton Board','Color','w',...
                 'Position',[200 100 700 850]);
    ax = axes('Parent',fig);
    hold(ax,'on');
    xlim(ax,[boardLeft-1, boardRight+1]);
    ylim(ax,[boardBottom-1, boardTop+2]);
    set(ax,'Color',[0.15 0.15 0.2],'XTick',[],'YTick',[],'Box','on');
    title(ax,'Galton Board Simulation','Color','w','FontSize',14);

    %% Draw pegs
    for level = 1:nLevels
        nPegs = level;
        y = boardTop - level * levelHeight;
        xOffset = -(nPegs - 1) / 2 * pegSpacing;
        for p = 1:nPegs
            x = xOffset + (p-1) * pegSpacing;
            plot(ax, x, y, '.', 'Color', [0.6 0.6 0.6], 'MarkerSize', 8);
        end
    end

    %% Draw bin walls
    binFloor = boardBottom;
    binTop   = boardBottom + 6;  % max visual height of bins
    for k = 0:nBins
        xWall = boardLeft + k * binWidth;
        plot(ax, [xWall xWall], [binFloor binTop], '-', 'Color', [0.4 0.4 0.4], 'LineWidth', 1);
    end
    plot(ax, [boardLeft, boardRight], [binFloor binFloor], '-', 'Color', [0.4 0.4 0.4], 'LineWidth', 1);

    %% Prepare bin counts and ball graphics in bins
    binCounts = zeros(1, nBins);
    binCenters = boardLeft + (0:nBins-1) * binWidth + binWidth/2;

    % Preallocate theoretical curve (drawn later)
    theorLine = plot(ax, NaN, NaN, 'r-', 'LineWidth', 2);

    % Counter text
    counterText = text(ax, boardRight - 1, boardTop + 1, '0 / 300', ...
                       'Color','w','FontSize',12,'HorizontalAlignment','right');

    %% Drop balls one by one
    ballRadius = 0.15;
    stepsPerLevel = 6;  % animation frames between peg levels

    for b = 1:nBalls
        % Determine path: at each level, go left (0) or right (1)
        choices = randi([0 1], 1, nLevels);
        binIndex = sum(choices) + 1;  % which bin (1..nBins)

        % Animate the ball falling through pegs
        xPos = 0;
        yPos = boardTop + 1;
        ballGfx = plot(ax, xPos, yPos, 'o', ...
                       'MarkerSize', 7, 'MarkerFaceColor', [0.2 0.6 1], ...
                       'MarkerEdgeColor', 'w', 'LineWidth', 0.5);

        for level = 1:nLevels
            % Target position after this peg
            xTarget = -(level-1)/2 * pegSpacing + (sum(choices(1:level))-0.5*choices(level)) * pegSpacing;
            if choices(level) == 1
                xTarget = xTarget + pegSpacing * 0.5;
            else
                xTarget = xTarget - pegSpacing * 0.5;
            end
            % Simpler: recompute from scratch
            xTarget = -(level)/2 * pegSpacing + sum(choices(1:level)) * pegSpacing;
            yTarget = boardTop - level * levelHeight - 0.3;

            xStart = xPos; yStart = yPos;
            for s = 1:stepsPerLevel
                t = s / stepsPerLevel;
                xPos = xStart + (xTarget - xStart) * t;
                yPos = yStart + (yTarget - yStart) * t;
                set(ballGfx, 'XData', xPos, 'YData', yPos);
                drawnow limitrate;
            end
        end

        % Animate ball falling into its bin
        xFinal = binCenters(binIndex);
        binCounts(binIndex) = binCounts(binIndex) + 1;
        yFinal = binFloor + (binCounts(binIndex) - 0.5) * (binWidth * 0.7);

        xStart = xPos; yStart = yPos;
        fallSteps = 8;
        for s = 1:fallSteps
            t = s / fallSteps;
            xPos = xStart + (xFinal - xStart) * t;
            yPos = yStart + (yFinal - yStart) * t;
            set(ballGfx, 'XData', xPos, 'YData', yPos);
            drawnow limitrate;
        end

        % Leave the ball in the bin (change color to stack nicely)
        set(ballGfx, 'XData', xFinal, 'YData', yFinal, ...
            'MarkerFaceColor', [0.3 0.75 1], 'MarkerSize', 6);

        % Update counter
        set(counterText, 'String', sprintf('%d / %d', b, nBalls));

        % Update theoretical normal curve overlay every 10 balls
        if mod(b, 10) == 0 || b == nBalls
            mu = nLevels / 2;
            sigma = sqrt(nLevels / 4);
            xCurve = linspace(0, nLevels, 200);
            yCurve = b * binWidth * normpdf(xCurve, mu, sigma);
            % Map to board coordinates
            xPlot = boardLeft + (xCurve + 0.5) * binWidth;
            yPlot = binFloor + yCurve * (binWidth * 0.7);
            set(theorLine, 'XData', xPlot, 'YData', yPlot);
        end

        % Adjust y-axis to fit tallest bin stack and Gaussian curve
        stackTop = max(binCounts) * binWidth * 0.7 + 1;
        % Also account for Gaussian curve peak height
        if b > 0
            mu = nLevels / 2;
            sigma = sqrt(nLevels / 4);
            curveTop = b * binWidth * normpdf(mu, mu, sigma) * (binWidth * 0.7);
        else
            curveTop = 0;
        end
        contentTop = max(stackTop, curveTop) + 1;
        yTop = max(boardTop + 2, binFloor + contentTop);
        ylim(ax, [boardBottom - 1, yTop]);
    end

    % Final annotation
    text(ax, 0, boardTop + 1.5, sprintf('Normal Distribution (n=%d, levels=%d)', nBalls, nLevels), ...
         'Color', 'y', 'FontSize', 11, 'HorizontalAlignment', 'center');

    hold(ax, 'off');
end
