This is working code written in Visual FoxPro, and object oriented version of xBase.
Return to A Study in Sphere Packing

*
*  packTurtle is a subclass of Turtle, previously
*  defined in procedure file 4DLOGO.PRG
*
*  Calculate IVM positions using quadrays and
*  output a PovRay (POV) file that will render
*  unit-radius spheres at these positions

*  By Kirby Urner, Oct 29, 1997
*  Last upgrade:

set procedure to 4dlogo  && make sure Turtle class is accessible

opacker = createobject("packTurtle")
opacker.makepoint(0,0,0,0)
opacker.calcve(1,0,0,0,0)
opacker.writepovray("ivm.pov",1)

release objects
set procedure to
return

define class packturtle as turtle
    dimension vevector(12,2)

    procedure init
        use datapoints order coords
        set safety off
        zap
        set safety on
        this.defineve()
        turtle::init()
    endproc

    procedure destroy
        use
    endproc

    procedure defineve
        local i,j,k
        k=1
        for i = 1 to 4
            for j = 1 to 4
                if i=j
                    loop
                endif
                this.vevector(k,1)=i
                this.vevector(k,2)=j
                k=k+1
            endfor
        endfor
    endproc

    procedure calcve(n,a,b,c,d)
        local k, q(4)
        for k = 1 to 12
            q=1
            q(this.vevector(k,1))=2
            q(this.vevector(k,2))=0
            q(1)=a+q(1)
            q(2)=b+q(2)
            q(3)=c+q(3)
            q(4)=d+q(4)
            this.makepoint(q(1),q(2),q(3),q(4))
            if n>1
                this.calcve(n-1,q(1),q(2),q(3),q(4))
            endif
        endfor
    endproc

    procedure makepoint(a,b,c,d)
        this.quadmov(a,b,c,d)
        this.storepos()
        this.quadmov(-a,-b,-c,-d)
    endproc

    procedure storepos(a,b,c,d)
        srch=str(int(this.quadpos(1)),2)+;
            str(int(this.quadpos(2)),2)+;
            str(int(this.quadpos(3)),2)+;
            str(int(this.quadpos(4)),2)
        if not seek(srch,"datapoints")
            insert into datapoints ;
                (acoord, bcoord, ccoord, dcoord, ;
                xcoord, ycoord, zcoord) ;
                values (this.quadpos(1), this.quadpos(2), ;
                this.quadpos(3), this.quadpos(4),;
                this.xyzpos(1), this.xyzpos(2), this.xyzpos(3))
        endif
    endproc

    procedure writepovray(filename,n)
        local hnd

        if file(filename)
            erase (filename)
        endif

        hnd=fcreate(filename)

        if hnd>0
            =fopen(filename)
        endif

        =fputs(hnd, "//POV-Ray script")
        =fputs(hnd, '#include"colors.inc"')
        =fputs(hnd, "")
        =fputs(hnd, "#declare Cam_factor = " + str(10 * n,3))
        =fputs(hnd, "#declare Camera_X = 1 * Cam_factor")
        =fputs(hnd, "#declare Camera_Y = 0.5 * Cam_factor")
        =fputs(hnd, "#declare Camera_Z = -0.7 * Cam_factor")
        =fputs(hnd, "camera { location  ")
        =fputs(hnd, "		up        <0, 1.0,  0>    right     <-1.33, 0,  0>")
        =fputs(hnd, "		direction <0, 0,  3>      look_at   <0, 0, 0> }")
        =fputs(hnd, "")
        =fputs(hnd, "light_source {  color White }")
        =fputs(hnd, "light_source {  color White }")
        =fputs(hnd, "")
        =fputs(hnd, "// Background:")
        =fputs(hnd, "background {color White}")

        =fputs(hnd, "#declare IVM=")
        =fputs(hnd, "union {")
        go top
        nocolorize = .t.
        scan while not eof()
            if nocolorize
                =fputs(hnd, "sphere{<"+str(xcoord,10,7)+",";
                    +str(ycoord,10,7)+",";
                    +str(zcoord,10,7)+">,1.00";
                    +" pigment {color Green }}")
            else
                =fputs(hnd, "sphere{<"+str(xcoord,10,7)+",";
                    +str(ycoord,10,7)+",";
                    +str(zcoord,10,7)+">,1.00";
                    +" pigment {color rgb<";
                    +str(rcolor,4,2)+",";
                    +str(gcolor,4,2)+",";
                    +str(bcolor,4,2)+">}}")
            endif
        endscan
        =fputs(hnd, "rotate <0, 0, 0> }")
        =fputs(hnd, "object{IVM}")
        =fclose(hnd)
    endproc

enddefine