POINT Software Programming/PPROG1.ISO/misc/euphoria/mset.ex

Text File | 1994-04-19 | 9.7 KB | 380 lines

        ---------------------------------
        -- Plotting the Mandelbrot Set --
        ---------------------------------
-- Usage: ex mset [filename.dat]
-- e.g. ex mset msetb

-- Either generate the initial picture, or redisplay an old one.
-- Hit Enter at any time to stop the display and then hit Enter again
-- to display a grid. Use the arrow keys to select the most interesting
-- box in the grid. Hit Enter to enlarge this box to the full size of
-- the screen, or hit q to quit. The pictures that you display are saved in
-- mseta.dat, msetb.dat, ... You can redisplay them and then put up the
-- grid and continue zooming in. As you zoom in, black areas are eroded
-- around the edges, as the iteration count is increased, and we find that
-- points originally believed to be members of the Mset are not members
-- after all.

without type_check

include graphics.e
include select.e
include get.e

constant GRAPHICS_MODE = 261 -- SVGA
                 -- choose a good mode for your machine
                 -- see euphoria\include\graphics.e

constant ZOOM_FACTOR = 20    -- grid size for zooming in

constant FALSE = 0, TRUE = 1
constant REAL = 1, IMAG = 2

constant ARROW_LEFT  = 331,
     ARROW_RIGHT = 333,
     ARROW_UP    = 328,
     ARROW_DOWN  = 336

    -- types --

type natural(integer x)
    return x >= 0
end type

type complex(sequence x)
    return length(x) = 2 and atom(x[1]) and atom(x[2])
end type

procedure beep()
-- make a beep sound
    atom t

    t = time()
    sound(500)
    while time() < t + .2 do
    end while
    sound(0)
end procedure

natural ncolors
integer max_color, min_color

object prev_mixture
prev_mixture = {0, 0, 0} -- color 0: black

procedure randomize_palette()
-- choose random color mixtures
    for i = max_color to min_color by -1 do
    prev_mixture = palette(i, rand(repeat(63, 3)))
    end for
    prev_mixture = palette(0, rand(repeat(63, 3)))
end procedure

procedure rotate_palette()
-- rotate the color mixtures of all the colors

    for i = max_color to min_color by -1 do
    prev_mixture = palette(i, prev_mixture)
    if atom(prev_mixture) then
        return -- didn't work
    end if
    end for
    prev_mixture = palette(0, prev_mixture)
end procedure

natural max_iter

sequence vc -- current video configuration

procedure grid(sequence x, sequence y, natural color)
-- draw the grid
    atom dx, dy

    dx = vc[VC_XPIXELS]/ZOOM_FACTOR
    dy = vc[VC_YPIXELS]/ZOOM_FACTOR

    for i = x[1] to x[2] do
    draw_line(color, {{i*dx, y[1]*dy}, {i*dx, y[2]*dy}})
    end for
    for i = y[1] to y[2] do
    draw_line(color, {{x[1]*dx, i*dy}, {x[2]*dx, i*dy}})
    end for
end procedure

function zoom()
-- select place to zoom in on next time
    integer key
    sequence box

    grid({0, ZOOM_FACTOR}, {0, ZOOM_FACTOR}, 7)
    box = {0, ZOOM_FACTOR-1}
    while TRUE do
    grid({box[1], box[1]+1}, {box[2], box[2]+1}, 14)
    key = get_key()
    if key != -1 then
        grid({box[1], box[1]+1}, {box[2], box[2]+1}, 7)
        if key = ARROW_UP then
        if box[2] > 0 then
            box[2] = box[2]-1
        end if
        elsif key = ARROW_DOWN then
        if box[2] < ZOOM_FACTOR-1 then
            box[2] = box[2]+1
        end if
        elsif key = ARROW_RIGHT then
        if box[1] < ZOOM_FACTOR-1 then
            box[1] = box[1]+1
        end if
        elsif key = ARROW_LEFT then
        if box[1] > 0 then
            box[1] = box[1]-1
        end if
        elsif key >= ' '  then
        return {}  -- quit
        else
        return {box[1], ZOOM_FACTOR - 1  - box[2]}
        end if
    end if
    end while
end function

procedure save_points(integer file, integer rep_count, integer color)
-- We do a bit of image compression here by recording the number of
-- consecutive points of a certain color. Can have up to 256 colors.
    while rep_count > 255 do
    puts(file, 255)
    puts(file, color)
    rep_count = rep_count - 255
    end while
    if rep_count > 0 then
    puts(file, rep_count)
    puts(file, color)
    end if
end procedure

function mset(complex lower_left,  -- lower left corner
          complex upper_right) -- upper right corner
-- Plot the Mandelbrot set over some region.
-- The Mandelbrot set is defined by the equation: z = z * z + C
-- where z and C are complex numbers. The starting point for z is 0.
-- If, for a given value of C, z approaches infinity, C is considered to
-- *not* be a member of the set. It can be shown that if the absolute value
-- of z ever becomes greater than 2, then the value of z will increase
-- towards infinity from then on. After a large number of iterations, if
-- the absolute value of z is still less than 2 then we assume with high
-- probability that C is a member of the Mset and this program will show
-- that point in black.
    complex c
    atom zr, zi, zr2, zi2, cr, ci, xsize, ysize
    natural member, stop, color, rep_count, width, height
    natural file_no
    integer pic_save, prev_color

    clear_screen()
    height = vc[VC_YPIXELS]
    width = vc[VC_XPIXELS]
    ncolors = vc[VC_NCOLORS]
    xsize = (upper_right[REAL] - lower_left[REAL])/(width - 1)
    ysize = (upper_right[IMAG] - lower_left[IMAG])/(height - 1)
    c = {0, 0}

    -- choose a new file to save the picture into
    file_no = 0
    for i = 'a' to 'z' do
    pic_save = open("mset" & i & ".dat", "rb")
    if pic_save = -1 then
        file_no = i
        exit
    else
        -- file exists
        close(pic_save)
    end if
    end for
    if file_no then
    pic_save = open("mset" & file_no & ".dat", "wb")
    else
    puts(1, "Couldn't find a new file name to use\n")
    return 1
    end if

    -- save graphics mode and max_iter
    printf(pic_save, "%d ", vc[VC_MODE])
    printf(pic_save, "%d ", max_iter)
    -- save lower_left & upper_right as floating-point sequences
    printf(pic_save, "{%20.15f,%20.15f}", lower_left)
    printf(pic_save, "{%20.15f,%20.15f}", upper_right)
    max_color = -1
    min_color = 99999
    for y = 0 to height - 1 do
    if get_key() != -1 then
        close(pic_save)
        return 0
    end if
    c[IMAG] = upper_right[IMAG] - y * ysize
    prev_color = -1 -- start fresh for each line
    rep_count = 0
    for x = 0 to width - 1 do
        c[REAL] = lower_left[REAL] + x * xsize
        member = TRUE
        zr = 0
        zi = 0
        zr2 = 0
        zi2 = 0
        cr = c[REAL]
        ci = c[IMAG]
        for i = 1 to max_iter do
        zi = 2 * zr * zi + ci
        zr = zr2 - zi2 + cr
        zr2 = zr * zr
        zi2 = zi * zi
        if zr2 + zi2 > 4 then
            member = FALSE
            stop = i
            exit
        end if
        end for
        if member = TRUE then
        color = 0
        else
        color = stop + 51 -- gives nice sequence of colors
        while color >= ncolors do
            color = color - ncolors
        end while
        if color > max_color then
            max_color = color
        end if
        if color < min_color then
            min_color = color
        end if
        end if
        pixel(color, {x, y}) -- also show non-member "bands"
        if color = prev_color then
        rep_count = rep_count + 1
        else
        save_points(pic_save, rep_count, prev_color)
        rep_count = 1
        prev_color = color
        end if
    end for
    -- close off count at end of each line
    save_points(pic_save, rep_count, color)
    end for
    beep()
    close(pic_save)
    return 0
end function

procedure view(integer pic_save)
-- redisplay a saved picture file
    integer x, color, rep_count

    max_color = -1
    min_color = 99999
    for y = 0 to vc[VC_YPIXELS] - 1 do
    x = 0
    while x < vc[VC_XPIXELS] do
        rep_count = getc(pic_save)
        color = getc(pic_save)
        if rep_count <= 0 then
        return
        end if
        if rep_count = 1 then
        pixel(color, {x, y})  -- faster
        x = x + 1
        else
        draw_line(color, {{x, y}, {x+rep_count-1, y}})
        x = x + rep_count
        end if
        if color != 0 then
        if color > max_color then
            max_color = color
        end if
        if color < min_color then
            min_color = color
        end if
        end if
    end while
    end for
end procedure

procedure Mandelbrot()
-- main procedure
    sequence delta, new_box
    complex lower_left, upper_right
    sequence cl, dataname, g
    integer pic_save, mode

    cl = command_line()
    if length(cl) >= 3 then
    -- redisplay a saved picture
    dataname = cl[3]
    pic_save = open(dataname, "rb")
    if pic_save = -1 then
        if not find('.', dataname) then
        dataname = dataname & ".dat"
        pic_save = open(dataname, "rb")
        end if
        if pic_save = -1 then
        puts(1, "Couldn't open " & dataname & '\n')
        return
        end if
    end if
    g = {}
    for i = 1 to 4 do
        g = g & get(pic_save)
    end for
    if g[1] != GET_SUCCESS or
       g[3] != GET_SUCCESS or
       g[5] != GET_SUCCESS or
       g[7] != GET_SUCCESS then
        puts(1, "Couldn't read " & dataname & '\n')
        return
    end if
    mode = g[2]
    max_iter = g[4]
    lower_left = g[6]
    upper_right = g[8]
    if graphics_mode(mode) then
    end if
    vc = video_config()
    view(pic_save)
    else
    -- initially show the upper half:
    max_iter = 30 -- increases as we zoom in
    lower_left = {-1, 0}
    upper_right = {1, 1}
    -- set up for desired graphics mode
    if not select_mode(GRAPHICS_MODE) then
        puts(2, "couldn't find a good graphics mode\n")
        return
    end if
    vc = video_config()
    if mset(lower_left, upper_right) then
        return
    end if
    end if

    while TRUE do
    while get_key() = -1 do
            --rotate_palette()
        randomize_palette()
    end while
    new_box = zoom()
    if length(new_box) = 0 then
        exit
    end if
    delta = (upper_right - lower_left)
    lower_left = lower_left + new_box / ZOOM_FACTOR * delta
    upper_right = lower_left + delta / ZOOM_FACTOR
    max_iter = max_iter * 2  -- need more iterations as we zoom in
    if mset(lower_left,  upper_right) then
        exit
    end if
    end while
end procedure

Mandelbrot()

if graphics_mode(-1) then
end if