home *** CD-ROM | disk | FTP | other *** search
/ Media Share 9 / MEDIASHARE_09.ISO / progmisc / euphor10.zip / MSET.EX < prev    next >
Text File  |  1993-06-18  |  10KB  |  381 lines

  1.         ---------------------------------
  2.         -- Plotting the Mandelbrot Set --
  3.         ---------------------------------
  4. -- Usage: ex mset [filename.dat]
  5. -- e.g. ex mset msetb
  6. -- Either generate the initial picture, or redisplay an old one.
  7. -- Hit Enter at any time to stop the display and then
  8. -- hit Enter again to display a grid.
  9. -- Use the arrow keys to select the most interesting box in the grid.
  10. -- Hit Enter to enlarge this box to the full size of the screen,
  11. -- or hit q to quit. The pictures that you display are saved in
  12. -- mseta.dat, msetb.dat, ...
  13. -- You can redisplay them and then put up the grid and continue zooming in.
  14. -- As you zoom in, black areas are eroded around the edges, as the iteration
  15. -- count is increased, and we find that points originally believed to
  16. -- be members of the Mset are not members after all.
  17.  
  18. without type_check
  19.  
  20. include graphics.e
  21. include select.e
  22. include get.e
  23.  
  24. constant GRAPHICS_MODE = 261 -- choose a good mode for your machine
  25.                  -- see euphoria\include\graphics.e
  26.  
  27. constant ZOOM_FACTOR = 20    -- grid size for zooming in
  28.  
  29. constant FALSE = 0, TRUE = 1
  30. constant REAL = 1, IMAG = 2
  31.  
  32. constant ARROW_LEFT  = 331,
  33.      ARROW_RIGHT = 333,
  34.      ARROW_UP    = 328,
  35.      ARROW_DOWN  = 336
  36.  
  37.     -- types --
  38.  
  39. type natural(integer x)
  40.     return x >= 0
  41. end type
  42.  
  43. type complex(sequence x)
  44.     return length(x) = 2 and atom(x[1]) and atom(x[2])
  45. end type
  46.  
  47. procedure beep()
  48. -- make a beep sound
  49.     atom t
  50.  
  51.     t = time()
  52.     sound(500)
  53.     while time() < t + .2 do
  54.     end while
  55.     sound(0)
  56. end procedure
  57.  
  58. natural ncolors
  59. integer max_color, min_color
  60.  
  61. object prev_mixture 
  62. prev_mixture = {0, 0, 0} -- color 0: black
  63.  
  64. procedure randomize_palette()
  65. -- choose random color mixtures    
  66.     for i = max_color to min_color by -1 do
  67.     prev_mixture = palette(i, rand(repeat(63, 3)))
  68.     end for
  69.     prev_mixture = palette(0, rand(repeat(63, 3)))
  70. end procedure
  71.  
  72. procedure rotate_palette()
  73. -- rotate the color mixtures of all the colors
  74.     
  75.     for i = max_color to min_color by -1 do
  76.     prev_mixture = palette(i, prev_mixture) 
  77.     if atom(prev_mixture) then
  78.         return -- didn't work
  79.     end if
  80.     end for
  81.     prev_mixture = palette(0, prev_mixture)
  82.     
  83. end procedure
  84.  
  85. natural max_iter
  86.  
  87. sequence vc -- current video configuration
  88.  
  89. procedure grid(sequence x, sequence y, natural color, natural thickness)
  90.     atom dx, dy
  91.  
  92.     dx = vc[VC_XPIXELS]/ZOOM_FACTOR
  93.     dy = vc[VC_YPIXELS]/ZOOM_FACTOR
  94.  
  95.     for i = x[1] to x[2] do
  96.     draw_line(color, thickness, {{i*dx, y[1]*dy}, {i*dx, y[2]*dy}})
  97.     end for
  98.     for i = y[1] to y[2] do
  99.     draw_line(color, thickness, {{x[1]*dx, i*dy}, {x[2]*dx, i*dy}})
  100.     end for
  101. end procedure
  102.  
  103. function zoom()
  104. -- draw a grid so we can pick a place to zoom in on next time
  105.     integer key
  106.     sequence box
  107.  
  108.     grid({0, ZOOM_FACTOR}, {0, ZOOM_FACTOR}, 7, 3)
  109.     box = {0, ZOOM_FACTOR-1}
  110.     while TRUE do
  111.     grid({box[1], box[1]+1}, {box[2], box[2]+1}, 14, 3)
  112.     key = get_key()
  113.     if key != -1 then
  114.         grid({box[1], box[1]+1}, {box[2], box[2]+1}, 7, 3)
  115.         if key = ARROW_UP then
  116.         if box[2] > 0 then
  117.             box[2] = box[2]-1
  118.         end if
  119.         elsif key = ARROW_DOWN then
  120.         if box[2] < ZOOM_FACTOR-1 then
  121.             box[2] = box[2]+1
  122.         end if
  123.         elsif key = ARROW_RIGHT then
  124.         if box[1] < ZOOM_FACTOR-1 then
  125.             box[1] = box[1]+1
  126.         end if
  127.         elsif key = ARROW_LEFT then
  128.         if box[1] > 0 then
  129.             box[1] = box[1]-1
  130.         end if
  131.         elsif key = 'q' then
  132.         return {}  -- quit
  133.         else
  134.         return {box[1], ZOOM_FACTOR - 1  - box[2]}
  135.         end if
  136.     end if
  137.     end while
  138. end function
  139.  
  140. procedure save_points(integer file, integer rep_count, integer color)
  141. -- record a number of points of a certain color in the file
  142. -- can only have up to 256 colors
  143.     while rep_count > 255 do
  144.     puts(file, 255)
  145.     puts(file, color)
  146.     rep_count = rep_count - 255
  147.     end while
  148.     if rep_count > 0 then
  149.     puts(file, rep_count)
  150.     puts(file, color)
  151.     end if
  152. end procedure
  153.  
  154. function mset(complex lower_left,  -- lower left corner
  155.           complex upper_right) -- upper right corner
  156. -- Plot the Mandelbrot set over some region.
  157. -- The Mandelbrot set is defined by the equation: z = z * z + C
  158. -- where z and C are complex numbers. The starting point for z is 0.
  159. -- If for a given value of C, z approaches infinity, C is considered to
  160. -- *not* be a member of the set. It can be shown that if the absolute value
  161. -- of z every becomes greater than 2, then the value of z will increase
  162. -- towards infinity from then on. After a large number of iterations, if
  163. -- the absolute value of z is still less than 2 then we assume with high
  164. -- probability that C is a member of the Mset and this program will show
  165. -- that point in black.
  166.     complex c
  167.     atom zr, zi, zr2, zi2, cr, ci, xsize, ysize
  168.     natural member, stop, color, rep_count, width, height
  169.     natural file_no
  170.     integer pic_save, prev_color
  171.  
  172.     clear_screen()
  173.     height = vc[VC_YPIXELS]
  174.     width = vc[VC_XPIXELS]
  175.     ncolors = vc[VC_NCOLORS]
  176.     xsize = (upper_right[REAL] - lower_left[REAL])/(width - 1)
  177.     ysize = (upper_right[IMAG] - lower_left[IMAG])/(height - 1)
  178.     c = {0, 0}
  179.  
  180.     -- choose a new file to save the picture into
  181.     file_no = 0
  182.     for i = 'a' to 'z' do
  183.     pic_save = open("mset" & i & ".dat", "rb")
  184.     if pic_save = -1 then
  185.         file_no = i
  186.         exit
  187.     else
  188.         -- file exists
  189.         close(pic_save)
  190.     end if
  191.     end for
  192.     if file_no then
  193.     pic_save = open("mset" & file_no & ".dat", "wb")
  194.     else
  195.     puts(1, "Couldn't find a new file name to use\n")
  196.     return 1
  197.     end if
  198.  
  199.     -- save graphics mode and max_iter
  200.     printf(pic_save, "%d ", vc[VC_MODE])
  201.     printf(pic_save, "%d ", max_iter)
  202.     -- save lower_left & upper_right as floating-point sequences
  203.     printf(pic_save, "{%20.15f,%20.15f}", lower_left)
  204.     printf(pic_save, "{%20.15f,%20.15f}", upper_right)
  205.     max_color = -1
  206.     min_color = 99999
  207.     for y = 0 to height - 1 do
  208.     if get_key() != -1 then
  209.         close(pic_save)
  210.         return 0
  211.     end if
  212.     c[IMAG] = upper_right[IMAG] - y * ysize
  213.     prev_color = -1 -- start fresh for each line
  214.     rep_count = 0
  215.     for x = 0 to width - 1 do
  216.         c[REAL] = lower_left[REAL] + x * xsize
  217.         member = TRUE
  218.         zr = 0
  219.         zi = 0
  220.         zr2 = 0
  221.         zi2 = 0
  222.         cr = c[REAL]
  223.         ci = c[IMAG]
  224.         for i = 1 to max_iter do
  225.         zi = 2 * zr * zi + ci
  226.         zr = zr2 - zi2 + cr
  227.         zr2 = zr * zr
  228.         zi2 = zi * zi
  229.         if zr2 + zi2 > 4 then
  230.             member = FALSE
  231.             stop = i
  232.             exit
  233.         end if
  234.         end for
  235.         if member = TRUE then
  236.         color = 0
  237.         else
  238.         color = stop + 51 -- gives nice sequence of colors
  239.         while color >= ncolors do
  240.             color = color - ncolors
  241.         end while
  242.         if color > max_color then
  243.             max_color = color
  244.         end if
  245.         if color < min_color then
  246.             min_color = color
  247.         end if
  248.         end if
  249.         pixel(color, {x, y}) -- also show non-member "bands"
  250.         if color = prev_color then
  251.         rep_count = rep_count + 1
  252.         else
  253.         save_points(pic_save, rep_count, prev_color)
  254.         rep_count = 1
  255.         prev_color = color
  256.         end if
  257.     end for
  258.     -- close off count at end of each line
  259.     save_points(pic_save, rep_count, color)
  260.     end for
  261.     beep()
  262.     close(pic_save)
  263.     return 0
  264. end function
  265.  
  266. procedure view(integer pic_save)
  267. -- redisplay a saved picture file
  268.     integer x, color, rep_count
  269.  
  270.     max_color = -1
  271.     min_color = 99999
  272.     for y = 0 to vc[VC_YPIXELS] - 1 do
  273.     x = 0
  274.     while x < vc[VC_XPIXELS] do
  275.         rep_count = getc(pic_save)
  276.         color = getc(pic_save)
  277.         if rep_count <= 0 then
  278.         return
  279.         end if
  280.         if rep_count = 1 then
  281.         pixel(color, {x, y})  -- faster
  282.         x = x + 1
  283.         else
  284.         draw_line(color, 1, {{x, y}, {x+rep_count-1, y}})
  285.         x = x + rep_count
  286.         end if
  287.         if color != 0 then
  288.         if color > max_color then
  289.             max_color = color
  290.         end if
  291.         if color < min_color then
  292.             min_color = color
  293.         end if
  294.         end if
  295.     end while
  296.     end for
  297. end procedure
  298.  
  299. procedure Mandelbrot()
  300. -- main procedure
  301.     sequence delta, new_box
  302.     complex lower_left, upper_right
  303.     sequence cl, dataname, g
  304.     integer pic_save, mode
  305.  
  306.     cl = command_line()
  307.     if length(cl) >= 3 then
  308.     -- redisplay a saved picture
  309.     dataname = cl[3]
  310.     pic_save = open(dataname, "rb")
  311.     if pic_save = -1 then
  312.         if not find('.', dataname) then
  313.         dataname = dataname & ".dat"
  314.         pic_save = open(dataname, "rb")
  315.         end if
  316.         if pic_save = -1 then
  317.         puts(1, "Couldn't open " & dataname & '\n')
  318.         return
  319.         end if
  320.     end if
  321.     g = {}
  322.     for i = 1 to 4 do
  323.         g = g & get(pic_save)
  324.     end for
  325.     if g[1] != GET_SUCCESS or
  326.        g[3] != GET_SUCCESS or
  327.        g[5] != GET_SUCCESS or
  328.        g[7] != GET_SUCCESS then
  329.         puts(1, "Couldn't read " & dataname & '\n')
  330.         return
  331.     end if
  332.     mode = g[2]
  333.     max_iter = g[4]
  334.     lower_left = g[6]
  335.     upper_right = g[8]
  336.     graphics_mode(mode)
  337.     vc = video_config()
  338.     view(pic_save)
  339.     else
  340.     -- initially show the upper half:
  341.     max_iter = 30 -- increases as we zoom in
  342.     lower_left = {-1, 0}
  343.     upper_right = {1, 1}
  344.     -- set up for desired graphics mode
  345.     if not select_mode(GRAPHICS_MODE) then
  346.         puts(2, "couldn't find a good graphics mode\n")
  347.         return
  348.     end if
  349.     vc = video_config()
  350.     if mset(lower_left, upper_right) then
  351.         return
  352.     end if
  353.     end if
  354.  
  355.     while TRUE do
  356.     while get_key() = -1 do
  357.             --rotate_palette()
  358.         randomize_palette()
  359.     end while
  360.     new_box = zoom()
  361.     if length(new_box) = 0 then
  362.         graphics_mode(3)
  363.         exit
  364.     end if
  365.     delta = (upper_right - lower_left)
  366.     lower_left = lower_left + new_box / ZOOM_FACTOR * delta
  367.     upper_right = lower_left + delta / ZOOM_FACTOR
  368.     max_iter = max_iter * 2  -- need more iterations as we zoom in
  369.     if mset(lower_left,  upper_right) then
  370.         exit
  371.     end if
  372.     end while
  373. end procedure
  374.  
  375. Mandelbrot()
  376. if vc[VC_COLOR] then
  377.     graphics_mode(3)
  378. else
  379.     graphics_mode(7)
  380. end if
  381.