Category Archives: programming

Lal, 1991 in situ 10-Be production rates

10Be is a cosmogenic radioactive nuclide that is produced when high energy cosmic rays collide with nuclides and cause spallation. 10Be is produced in the atmosphere (and then transported down to the surface) as “meteoric”, and produced within mineral lattices in soil and rocks as “in situ“. In 1991, Devendra Lal wrote a highly cited paper about the physics of in situ produced Beryllium-10 (10Be). In the paper he lays out an equation for the production of in situ 10Be (q) based on latitude and altitude. I’m currently working on an idea I have for using cosmogenic nuclides as tracers for basin scale changes in uplift rate, so I wanted to see what his equation looked like applied. The equation is a third degree polynomial, with coefficients that depend on latitude (L), and direct dependency on altitude (y).

I grabbed an old raster (GEBCO 2014 30 arc second) I had laying around for Eastern North America and plotted it up. First, the elevation map (obviously latitude is on the y-axis…)


Elevation map for ENAM.

And then apply the Lal, 1991 equation and find


Plotting Lal’s 1991 in situ production rate equation for ENAM. Green–>red increasing production rate. Production rate = NA in water.

I think the interesting observation is for how little of the mapped area there is any significant change in the production rate. Maybe this should be obvious since the polynomial has direct dependence on altitude and altitude doesn’t change that much in most of the map. Further the dependence of latitude is not all all observable with this map; perhaps because the latitude range is not very large, or the coefficients never change by more than an order of magnitude anyway. Next time, maybe a world elevation map! Not sure my computer has enough memory…

You can grab the code I used from here and Lal’s paper from here.

River hysteresis and plotting hysteresis data in R

Hysteresis is the concept that a system (be it mechanical, numerical, or whatever else) is dependent of the history of the system, and not only the present conditions. This is the case for rivers. For example, consider the following theoretical flood curve and accompanied sediment discharge curve (Fig. 1a).

Figure 1. Theoretical plots to demonstrate the hysteresis of a river.

Figure 1. Theoretical plots to demonstrate the hysteresis of a river.

With the onset of the flood, the increased sediment transport capacity of the system entrains more sediment and the sediment discharge curve (red, Fig. 1a) rises. However, the system may soon run out of sediment to transport (really just a reduction in easily transportable sediment), and the sediment discharge curve decreases although the water discharge curve remains high in flood.

In Fig. 1b, the sediment discharge and water discharge are plotted through time, a typical way of observing the hysteresis of a system. Note that for the rising limb and falling limb of the river flood, the same water discharge produces two different sediment transport values.

Now, let’s imagine that we want to investigate how important the history of the system is to the present state of our study river. You can grab the data I’ll use, here. This is data from one year of flow on the Huanghe (Yellow River) in China, and it has been smoothed with a moving average function to make the hysteresis function more visible.

Making the plot

It is easy enough to plot a line with R (the lines function) but with a hysteresis plot, it is important to be able to determine which direction time is moving forward along the curve. For this reason we want to use arrows. So we plot the line data first, with:


and then using a constructed vector of every 22nd number, we plot an arrow over top of the lines using:

s <- seq(from=1, to=length(df[,"Qw"])-1, by=22)
arrows(df[,"Qw"][s], df[,"Qs"][s], df[,"Qw"][s+1], df[,"Qs"][s+1], 

Finally, with a few more modifications to the plot (see the full code here), we can produce Fig. 2 below. This plot is comparable to the theoretical one above.

Figure 2: Hysteresis plot from actual data.

Figure 2. Hysteresis plot from actual data.

Using the green lines and points, I have highlighted the observation that for the rising limb and falling limb of a flood, there can be substantially different sediment discharges for the same water discharge — this observation is not so easily made from the plot on the left, but it is immediately clear in the hysteresis plot on the right.

Be inspired while coding! — Matlab script

Here is a joke script I wrote a while ago that returns an inspirational quote when called. I made it as a joke to send to my lab group during a particularly grueling week of coding. You can simply call the function at the beginning of your script to have a quote printed to the stdout or you can wrap it inside a wbar if you want.

Grab the code from my GitHub, here.


example output from the inspire.m script.

Quick script for connecting to University VPN

I can be pretty lazy at times, so much that I will go out of my way a bit to write a bit of code to simplify my life. I don’t often connect to my university’s VPN, but when I eventually do want to, I can never remember the command needed to do it. So I have to take the time and look it up. Well I cut that out today with a basic little bash script with a name I can remember when I need it. I titled my script “” but anything memorable to you would be fine. The script takes one argument (“on” or “off”) to either connect or disconnect with the proxy. Stick the script in your bin folder and it will execute simply with on.

The script is reproduced below, or you can grab it from here.

# Andrew J. Moodie
# Feb 2015
# easily remembered interface for vpn connection
if [ "$1" = "on" ]
    sudo /usr/sbin/vpnc RiceVPN.conf
elif [ "$1" = "off" ]
    echo " input error -- one argument must be given"
    echo " valid inputs: 'on' or 'off'"

Note that you could further automate this if you wanted to (e.g. auto-enter passwords, connect on startup), but this is a bad idea for both security and for internet speed.