Generating triangular/hexagonal coordinates (xyz)

I'm trying to come up with an iterative function that generates xyz coordinates for a hexagonal grid. With a starting hex position (say 0,0,0 for simplicity), I want to calculate the coordinates for each successive "ring" of hexagons, as illustrated here:

So far, all I've managed to come up with is this (example in javascript):

var radius = 3
var xyz = [0,0,0];

// for each ring
for (var i = 0; i < radius; i++) {
    var tpRing = i*6;
    var tpVect = tpRing/3;
    // for each vector of ring
    for (var j = 0; j < 3; j++) {
        // for each tile in vector
        for(var k = 0; k < tpVect; k++) {
            xyz[0] = ???;
            xyz[1] = ???;
            xyz[2] = ???;
            console.log(xyz);
        }
    }
}

I know each ring contains six more points than the previous and each 120° vector contains one additional point for each step from the center. I also know that x + y + z = 0 for all tiles. But how can I generate a list of coordinates that follow the sequence below?

    0, 0, 0

    0,-1, 1
    1,-1, 0
    1, 0,-1
    0, 1,-1
   -1, 1, 0
   -1, 0, 1

    0,-2, 2
    1,-2, 1
    2,-2, 0
    2,-1,-1
    2, 0,-2
    1, 1,-2
    0, 2,-2
   -1, 2,-1
   -2, 2, 0
   -2, 1, 1
   -2, 0, 2
   -1,-1, 2

Another possible solution, that runs in O(radius²), unlike the O(radius⁴) of tehMick's solution (at the expense of a lot of style) is this:

radius = 4
for r in range(radius):
    print "radius %d" % r
    x = 0
    y = -r
    z = +r
    print x,y,z
    for i in range(r):
        x = x+1
        z = z-1
        print x,y,z
    for i in range(r):
        y = y+1
        z = z-1
        print x,y,z
    for i in range(r):
        x = x-1
        y = y+1
        print x,y,z
    for i in range(r):
        x = x-1
        z = z+1
        print x,y,z
    for i in range(r):
        y = y-1
        z = z+1
        print x,y,z
    for i in range(r-1):
        x = x+1
        y = y-1
        print x,y,z

or written a little more concisely:

radius = 4
deltas = [[1,0,-1],[0,1,-1],[-1,1,0],[-1,0,1],[0,-1,1],[1,-1,0]]
for r in range(radius):
    print "radius %d" % r
    x = 0
    y = -r
    z = +r
    print x,y,z
    for j in range(6):
        if j==5:
            num_of_hexas_in_edge = r-1
        else:
            num_of_hexas_in_edge = r
        for i in range(num_of_hexas_in_edge):
            x = x+deltas[j][0]
            y = y+deltas[j][1]
            z = z+deltas[j][2]            
            print x,y,z

It's inspired by the fact the hexagons are actually on the exterior of a hexagon themselves, so you can find the coordinates of 1 of its points, and then calculate the others by moving on its 6 edges.