## Simulate a drop of liquid percolating down through a mixture of
## sand and rock.
##
## Generate a random 100x100 matrix of 0's and 1's where 1 corresponds
## to rock and 0 corresponds to sand.  Generate the values in the
## matrix independently, based on a given probability p that each
## space is occupied by rock.
##
## A drop of liquid starts in the middle of the top layer (row 1,
## column 50).  It then moves according to the following four options,
## where options with lower numbers have higher precedence.
##
## 1. If the space directly below is sand, move there.
## 2. If the space below and to the left is sand, move there.
## 3. If the space below and to the right is sand, move there.
## 4. If the space directly to the right is sand, move there.
##
## If none of these moves can be made, the drop of liquid is stuck.
##
## Use simulation to calculate the average depth to which the liquid
## drops before getting stuck, and the proportion of the time that the
## drop reaches the bottom layer.


## The density of rocks in the sand.
p = 0.5

## The number of simulation replications.
nrep = 1e3

## The total depth across the simulation replications.
TD = 0

## The number of times that the bottom is reached.
NB = 0

## Simulation replications.
for (k in 1:nrep) {

    ## Randomly lay out the rocks.
    M = 1*(runif(100*100) < p)
    M = matrix(M, nrow=100, ncol=100)

    ## The initial position of the droplet.
    r = 1
    c = 50

    ## Let the droplet percolate through the rocks.
    while (r<100) {

        ## Always go straight down if possible.
        if (M[r+1,c] == 0) {
            r = r+1
        }

        ## Next try down/left.
        else if ( (c>1) & (M[r+1,c-1] == 0) ) {
            r = r+1
            c = c-1
        }

        ## Next try down/right.
        else if ( (c<100) & (M[r+1,c+1] == 0) ) {
            r = r+1
            c = c+1
        }

        ## Next try right.
        else if ( (c<100) & (M[r,c+1] == 0) ) {
            c = c+1
        }

        ## We're stuck
        else { break } 
    }

    ## Keep track of how often we reach the bottom. 
    if (r==100) { NB = NB + 1 }

    ## Keep track of the total of the final depths.
    TD = TD + r
}

## The estimated probability that we reach the bottom.
NB = NB/nrep

## The average depth that is reached.
TD = TD/nrep