versioninfo()

# ジュリア集合を計算
function julia_set(c::Complex128)
    SIZE = 500
    x = linspace(-1.5, 1.5, SIZE)
    y = linspace(-1.0, 1.0, SIZE)
    z = zeros(Complex128, SIZE, SIZE)
    for j = 1:SIZE, i = 1:SIZE
        z[i, j] = x[j] + y[i] * im
    end
    k = 2.0
    julia = zeros(Float64, SIZE, SIZE)
    for i = 1:length(z), n = 1:500
        z[i] = z[i]^k + c
        if abs(z[i]) > 2.0
            julia[i] = 1.0 / (10.0+n)
            break
        end
    end
    julia
end

# 代表的なジュリア集合6つをプロット
function demo_julia_set()
    c = [Complex128(1 - (1+sqrt(5))/2),
          0.28500 + 0.0000im,
          0.28500 + 0.0100im,
         -0.70176 - 0.3842im,
         -0.83500 - 0.2321im,
         -0.80000 + 0.1560im,]

    clf()
    for i = 1:length(c)
        @time z = julia_set(c[i])
        @show minimum(z), maximum(z)
        subplot(2,3,i)
        show_contourf(z)
        title(string("c = ", c[i]))
        xticks([])
        yticks([])
    end
    subplots_adjust(left=0.05, bottom=0.05, right=0.95, top=0.95,
                    wspace=0.1, hspace=0.1)
end

function show_contourf(z)
    levels = matplotlib[:ticker][:MaxNLocator](nbins=50)[:tick_values](minimum(z), maximum(z))
    plt[:contourf]( z, levels=levels, cmap=get_cmap("gray_r"))
end

using PyPlot

demo_julia_set()

# i=j=0 # a1
for i=1:3, j=1:3
    @show i,j     
    if j==2 
        break
    end
end
# @show i, j # a2

    zc = Complex( 0.28500 , 0.0100 )
    z = [ Complex(x,y) for x in linspace(-0.5, 0.5, 3 ), y in linspace(-0.5, 0.5, 3 ) ]
    julia = zeros(Float64, size(z))
    nmax = 500
    for i = eachindex(z)
        julia[i] = nmax
        @time for n = 1:500
            if abs2(z[i]) > 4.0
                julia[i] = n
                break
            end
            z[i] = z[i] * z[i] + zc
        end
    end
    julia

    zc = Complex( 0.28500 , 0.0100 )
    const zs = [ Complex(x,y) for x in linspace(-0.5, 0.5, 3 ), y in linspace(-0.5, 0.5, 3 ) ]
    nmax = 500
    julia = zeros(Float64, size(zs))
    julia[:] = nmax
    for i = eachindex(zs)
        zz=zs[i]
        @time for n = 1:nmax
            if abs2(zz) > 4.0
                julia[i]=n
                break
            end
            zz = zz*zz + zc
        end
    end
    julia

function compute_julia_set(zc,ju,zs,nmax)
    @assert size(ju)==size(zs) "-- ju: size error"
    @inbounds ju[:] = 1.0 / (10.0+nmax)
    @time for i=eachindex(zs)
        zz=zs[i]
        for n=1:nmax
            if abs2(zz) > 4.0
                ju[i]=1.0 / (10.0+n)
                break
            end
            zz=zz*zz + zc
        end
    end
end

function do_compute_julia_set()
    const zs=[ Complex(x,y) for 
        x in linspace(-1.5, 1.5, 600 ), 
        y in linspace(-1.0, 1.0, 400 ) ]
    const xs=real(zs)
    const ys=imag(zs)
    @time ju=zeros(Float64, size(zs))
    fig = matplotlib[:pyplot][:gcf]()
    fig[:set_size_inches](18.5, 10.5)
    nmax=500
    result=[]
    nc=1
    for zc in [ 
            Complex( 1.0 - (1.0+sqrt(5.0))/2.0, 0.0),
            Complex( 0.28500,  0.0000 ),
            Complex( 0.28500,  0.0100 ),
            Complex(-0.70176, -0.3842 ),
            Complex(-0.83500, -0.2321 ),
            Complex(-0.80000,  0.1560 ) ]
        @time compute_julia_set( zc, ju, zs, nmax )
        # println( minimum(ju) )
        # println( maximum(ju) )
        push!( result, deepcopy(ju) )
        levels = matplotlib[:ticker][:MaxNLocator](nbins=50)[:tick_values](minimum(ju), maximum(ju))
        plt[:subplot](2,3,nc)
        plt[:contourf](xs,ys,ju, levels=levels, cmap=get_cmap("gray_r"))
        plt[:title](string("c = ", zc ))
        plt[:xticks]([])
        plt[:yticks]([])
        nc = nc + 1
    end
    # println( typeof(result) )
    # println( size(result) )
    plt[:subplots_adjust](left=0.05, bottom=0.05, right=0.95, top=0.95, wspace=0.1, hspace=0.1)
    # savefig("result6.jpg")
end

using PyPlot

do_compute_julia_set()