sandpile.java

模拟砂堆的元胞自动机源代码。对于认识和学习元胞自动机（CA）很有帮助

// Sandpile.java

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

public class Sandpile extends Animation {

    int L = 20;                 // length of pile
    int grains;                 // total grains added
    int lost;                   // grains lost off edge
    int[] height;               // height
    int[] slope;                // local slope
    int[] change;               // change in height in toppling event
    int criticalSlope = 1;      // topples if slope > critSlope
    boolean random;             // add grain at random position
    boolean unstable;           // pile is unstable
    boolean[] move;             // move unstable grains to neighbor sites
    int t;                      // time step
    int[] N;                    // accumulate number of topplings
    int[] E;                    // accumulate avalanche sizes
    int maxHeight;              // for plotting

    void initial () {
        if (height == null || height.length < L) {
            height = new int[L];
            slope = new int[L];
            change = new int[L];
            move = new boolean[L];
            N = new int[L];
            E = new int[L];
        }
        for (int i = 0; i < L; i++) {
            height[i] = slope[i] = change[i] = 0;
            move[i] = false;
        }
        for (int i = 0; i < N.length; i++)
            N[i] = 0;
        for (int i = 0; i < E.length; i++)
            E[i] = 0;
        maxHeight = L;
        unstable = false;
        grains = lost = 0;
        t = 0;
    }

    int topplings;
    int unstableSites;

    void checkStability () {
        unstable = false;
        for (int i = 0; i < L; i++) {
            if (slope[i] > criticalSlope) {
                unstable = move[i] = true;
                ++unstableSites;
            } else {
                move[i] = false;
            }
        }
        if (unstable) {
            ++topplings;
        } else {
            if (unstableSites > E.length - 1) {
                int len = E.length;
                int[] temp = new int[unstableSites + 10];
                for (int i = 0; i < len; i++)
                    temp[i] = E[i];
                E = temp;
            }
            ++E[unstableSites];
            unstableSites = 0;
            if (topplings > N.length - 1) {
                int len = N.length;
                int[] temp = new int[topplings + 10];
                for (int i = 0; i < len; i++)
                    temp[i] = N[i];
                N = temp;
            }
            ++N[topplings];
            topplings = 0;
            ++t;
        }
    }

    void timeStep () {
        for (int i = 0; i < L; i++) {
            height[i] += change[i];
            change[i] = 0;
            if (height[i] > maxHeight)
                maxHeight = height[i];
        }
        for (int i = 0; i < L - 1; i++)
            slope[i] = height[i] - height[i + 1];
        slope[L - 1] = height[L - 1];
        checkStability();
        if (unstable) {
            for (int i = 0; i < L; i++) {
                if (move[i]) {
                    change[i] -= criticalSlope;
                    unstableSites += criticalSlope;
                    for (int j = 1; j <= criticalSlope; j++) {
                        if (i + j < L)
                            change[i + j] += 1;
                        else
                            ++lost;
                    }
                }
            }
        } else {                // add grain
            int site = 0;
            if (random)
                site = (int) (L * Math.random());
            change[site] = 1;
            ++grains;
        }
    }

    class Pile extends Plot {
        Pile () {
            setSize(300, 300);
        }
        
        public void paint () {
            clear();
            setWindow(0, maxHeight + 1, 0, maxHeight + 1);
            setColor("yellow");
            for (int i = 0; i < L; i++) {
                for (int j = 0; j < outer.height[i]; j++) {
                    if (j == outer.height[i] + change[i])
                        break;
                    boxArea(i + 0.1, i + 0.9, j + 0.1, j + 0.9);
                }
            }
            setColor("red");
            for (int i = 0; i < L; i++) {
                for (int j = outer.height[i]; j < outer.height[i] + change[i]; j++)
                    boxArea(i + 0.1, i + 0.9, j + 0.1, j + 0.9);
            }
        }
    }

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

        // power law plot 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, 2.25);
            // topplings
            setColor("green");
            analyze(N);
            if (maxIndex == 0)
                return;
            for (int i = 0; i <= maxIndex; i++) {
                double x1 = i / (maxIndex + 1.0);
                double x2 = (i + 1) / (maxIndex + 1.0);
                double y1 = 1.1;
                double y2 = 1.1 + logData[i];
                boxArea(x1, x2, y1, y2);
            }
            plotStringLeft("power = " + Easy.format(exponent, 3), 1, 2.0);
            setColor("black");
            plotString("log(N)", 0, 2.15);
            plotLine(0, 1.1, 1, 1.1);
            plotLine(0, 1.1, 0, 2.1);
            double d = 0.5 / (maxIndex + 1);
            plotStringCenter("0", d, 1);
            plotStringCenter("Topplings", 0.5, 1);
            plotStringCenter("" + maxIndex, 1 - d, 1);
            // unstable sites
            setColor("blue");
            analyze(E);
            if (maxIndex == 0)
                return;
            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];
                boxArea(x1, x2, y1, y2);
            }
            plotStringLeft("power = " + Easy.format(exponent, 3), 1, 0.9);
            setColor("black");
            d = 0.5 / (maxIndex + 1);
            plotLine(0, 0, 1, 0);
            plotLine(0, 0, 0, 1);
            plotStringCenter("0", d, -0.1);
            plotStringCenter("Unstable sites", 0.5, -0.1);
            plotStringCenter("" + maxIndex, 1 - d, -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);
            plotStringCenter("Grains added " + grains, 30, 1);
            plotStringCenter("Grains lost " + lost, 50, 1);
            plotString("Topplings " + topplings, 65, 1);
            plotString("Unstable sites " + unstableSites, 80, 1);
        }
    }

    Pile pile;
    PowerLaw powerLaw;
    Output output;
    Reader LReader, criticalSlopeReader, skipReader;
    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(criticalSlopeReader = new Reader("Critical slope", criticalSlope));
        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();
            }
        });
        addControlPanel();
    }

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

    public void reset () {
        L = LReader.readInt();
        criticalSlope = criticalSlopeReader.readInt();
        skip = skipReader.readInt();
        initial();
        pile.repaint();
        powerLaw.repaint();
        output.repaint();
    }

    Sandpile outer;
    public Sandpile () {
        super();
        outer = this;
    }

    public static void main (String[] args) {
        Sandpile sandpile = new Sandpile();
        sandpile.frame("One-dimensional sandpile", 650, 500);
    }
}