Lecture 39: April 27

Two-dimensional sandpile automaton

Here is the code which generates this applet:

// Sand2D.java

import comphys.graphics.*;
import comphys.Easy;
import java.awt.*;
import java.awt.event.*;

public class Sand2D extends Animation {

    // 2-dimensional sandpile automaton

    int L = 21;                 // number of sites = L * L
    int[][] m;                  // slope at each site
    int[][] mOld;               // old value of slope at each site
    boolean[][] unstable;       // if m > 3
    boolean pileUnstable;       // at least one unstable site
    int t;                      // time step
    int topplings;              // number of unstable sites
    int[] N;                    // distribution of avalanche sizes
    int minSize;                // minimum size to determine exponent
    int maxSize;                // maximum size to determine exponent
    boolean random;             // if false add sand at central site

    void setBoundaryValues () {
        // open boundaries
        for (int i = 1; i <= L; i++)
            m[i][0] = m[0][i] = m[i][L + 1] = m[L + 1][i] = 0;
    }
    
    void initial () {
        if (m == null || m.length < L + 2) {
            m = new int[L + 2][L + 2];
            mOld = new int[L + 2][L + 2];
            unstable = new boolean[L + 2][L + 2];
            N = new int[L];
        }
        // interior values
        for (int i = 1; i <= L; i++) {
            for (int j = 1; j <= L; j++)
                m[i][j] = 0;
        }
        setBoundaryValues();
        for (int i = 0; i < L + 2; i++) {
            for (int j = 0; j < L + 2; j++)
                unstable[i][j] = false;
        }
        t = 0;
        topplings = 0;
        for (int i = 0; i < N.length; i++)
            N[i] = 0;
        minSize = 1;
        maxSize = L;
    }

    void addGrain () {
        int i = (L + 1) / 2;
        int j = (L + 1) / 2;
        if (random) {
            i = 1 + (int) (Math.random() * L);
            j = 1 + (int) (Math.random() * L);
        }
        m[i][j] += 1;
    }

    void markUnstableSites () {
        pileUnstable = false;
        for (int i = 1; i <= L; i++) {
            for (int j = 1; j <= L; j++) {
                if (m[i][j] > 3) {
                    unstable[i][j] = true;
                    pileUnstable = true;
                    ++topplings;
                } else {
                    unstable[i][j] = false;
                }
            }
        }
    }

    void topple () {
        for (int i = 1; i <= L; i++) {
            for (int j = 1; j <= L; j++) {
                if (unstable[i][j]) {
                    m[i][j] -= 4;
                    m[i][j - 1] += 1;
                    m[i][j + 1] += 1;
                    m[i - 1][j] += 1;
                    m[i + 1][j] += 1;
                }
            }
        }
    }

    void updatePile () {
        if (pileUnstable) {
            markUnstableSites();
            topple();
            setBoundaryValues();
        } else {
            if (N.length < topplings + 1) {
                int[] temp = new int[topplings + 10];
                for (int i = 0; i < N.length; i++)
                    temp[i] = N[i];
                N = temp;
            }
            ++N[topplings];
            topplings = 0;
            addGrain();
            ++t;
            markUnstableSites();
        }
    }

    class Pile extends Plot {
        Pile () {
            setSize(300, 300);
        }

        String[] color = {"white", "yellow", "orange", "red",
                          "green", "cyan", "blue", "magenta"};

        public void paint () {
            clear();
            setWindow(0, L + 2, 0, L + 2);
            for (int i = 0; i < L + 2; i++) {
                for (int j = 0; j < L + 2; j++) {
                    setColor(color[m[i][j]]);
                    boxArea(i, i + 1, j, j + 1);
                }
            }
        }
    }

    class PowerLaw extends Plot {
        PowerLaw () {
            setSize(300, 300);
        }

        // power law and exponent
        int maxIndex;
        double[] logData;
        double exponent;
        void analyze (int[] data) {
            if (logData == null || logData.length < data.length)
                logData = new double[data.length];
            maxIndex = 0;
            int numValues = 0;
            double maxValue = 0;
            double xAv = 0, yAv = 0, xxAv = 0, xyAv = 0;
            for (int i = 0; i < data.length; i++) {
                double value = 0;
                if (data[i] > 0) {
                    maxIndex = i;
                    value = Math.log(data[i]);
                    if (value > maxValue)
                        maxValue = value;
                    xAv += i;
                    yAv += value;
                    xxAv += i * (double) i;
                    xyAv += i * value;
                    ++numValues;
                } else {
                    logData[i] = 0;
                }
                logData[i] = value;
            }
            exponent = 0;
            if (numValues > 0) {
                for (int i = 0; i <= maxIndex; i++) {
                    logData[i] /= maxValue;
                }
                xAv /= numValues;
                yAv /= numValues;
                xxAv /= numValues;
                xyAv /= numValues;
                exponent = (xyAv - xAv * yAv) / (xxAv - xAv * xAv);
            }
        }

        public void paint () {
            clear();
            setWindow(-0.1, 1.1, -0.15, 1.15);
            setColor("cyan");
            analyze(N);
            double d = 0.5 / (maxIndex + 1);
            if (maxIndex > 0) {
                for (int i = 0; i <= maxIndex; i++) {
                    double x1 = i / (maxIndex + 1.0);
                    double x2 = (i + 1) / (maxIndex + 1.0);
                    double y1 = 0;
                    double y2 = logData[i];
                    plotLine(x2, y1, x2, y2);
                    boxArea(x1, x2, y1, y2);
                }
                setColor("blue");
                plotStringLeft("power = " + Easy.format(exponent, 3), 1, 1.0);
                plotStringCenter("0", d, -0.1);
                plotStringCenter("" + maxIndex, 1 - d, -0.1);
            }
            setColor("black");
            plotString("log(N)", 0, 1.05);
            plotLine(0, 0, 0, 1);
            plotLine(0, 0, 1, 0);
            plotStringCenter("Topplings", 0.5, -0.1);
        }
    }

    class Output extends Plot {
        Output () {
            setSize(610, 30);
        }

        public void paint () {
            clear();
            setWindow(0, 100, 0, 3);
            setColor("blue");
            boxArea(0, 100, 0, 3);
            setColor("white");
            plotString("time step = " + t, 5, 1);
            plotString("topplings = " + topplings, 25, 1);
        }
    }

    Pile pile;
    PowerLaw powerLaw;
    Output output;
    Reader LReader, skipReader, minSizeReader, maxSizeReader;
    Checkbox randomBox;
    int skip;

    public void init () {
        initial();
        add(pile = new Pile());
        add(powerLaw = new PowerLaw());
        add(output = new Output());
        add(LReader = new Reader("L = ", L));
        add(skipReader = new Reader("Skip steps", skip));
        add(randomBox = new Checkbox("Add random", random));
        randomBox.addItemListener(new ItemListener() {
            public void itemStateChanged(ItemEvent ie) {
                random = randomBox.getState();
            }
        });
        add(minSizeReader = new Reader("Min size = ", minSize));
        add(maxSizeReader = new Reader("Max size = ", maxSize));
        addControlPanel();
    }

    public void step () {
        for (int s = 0; s <= skip; s++)
            updatePile();
        pile.repaint();
        powerLaw.repaint();
        output.repaint();
    }

    public void reset () {
        L = LReader.readInt();
        skip = skipReader.readInt();
        minSize = minSizeReader.readInt();
        maxSize = maxSizeReader.readInt();
        initial();
        pile.repaint();
        powerLaw.repaint();
        output.repaint();
   }

    public static void main (String[] args) {
        Sand2D sand2D = new Sand2D();
        sand2D.frame("2-D Sandpile Automaton", 630, 600);
    }
}