mandelbulb.jl

using Plots, GLPlot; GLPlot.init()
using GLPlot, Reactive, GeometryTypes, Colors
using CUDAnative, GPUArrays
import GPUArrays: GPUArray, GLBackend, CUBackend
cuctx = CUBackend.init()

# Soon this should be fixed, but so far we need to qualify inbuild math functions with cu
const cu = CUDAnative

function mandelbulb{T}(x0::T,y0::T,z0::T, n, iter)
    x,y,z = x0,y0,z0
    for i=1:iter
        r = cu.sqrt(x*x + y*y + z*z) 
        theta = cu.atan2(cu.sqrt(x*x + y*y) , z)
        phi = cu.atan2(y,x)
        rn = cu.pow(r, n)
        x1 = rn * cu.sin(theta*n) * cu.cos(phi*n) + x0
        y1 = rn * cu.sin(theta*n) * cu.sin(phi*n) + y0
        z1 = rn * cu.cos(theta*n) + z0
        (x1*x1 + y1*y1 + z1*z1) > n && return T(i)
        x,y,z = x1,y1,z1
    end
    T(iter)
end

dims = (200,200,200)
vol = GPUArray(Float32, dims)
xrange, yrange, zrange = ntuple(3) do i
    # linearly spaced array (not dense) from -1 to 1
    # we need the reshape to make it work with broadcast
    d = dims[i]
    reshape(linspace(-1f0, 1f0, d), ntuple(j-> j==i ? d : 1, 3))
end
# create two sliders to interact with the 2 parameters of the mandelbulb function
itslider = GLPlot.play_widget(1:15);
nslider = GLPlot.play_widget(linspace(1f0,30f0, 100));
# register a callback to the sliders with foldp ("fold over past"), with v0 being the `past` array
volume = foldp(Array(vol), nslider, itslider) do v0, n, it
    # julia's new inplace broadcast syntax does wonders for GPU programming.
    # since vol is a GPUArray, this will dispatch to the GPU implementation of broadcast!
    # Since we can compile a large set of julia code for the GPU, `mandelbulb` can be broadcastet like this!
    vol .= mandelbulb.(xrange, yrange, zrange, n, it)
    unsafe_copy!(v0, vol) # copy to RAM. In future versions of GLPlot/GLVisualize, this won't be necessary
    v0
end

glplot(volume)