examples/mandelbrot.4th

   1 ( --- Complex arithmetic --- )
   2
   3 ( Location of floating point. )
   4 10000 value precision
   5
   6 : sign dup abs / ;
   7
   8 : >scaled
   9     precision 10 / * over
  10     ?dup 0<> if
  11         sign *
  12     then
  13     swap precision * +
  14 ;
  15
  16 ( Redefine multiplication.  Yay forth! )
  17 : * precision */ ;
  18
  19 : c* ( x1 y1 x2 y2 -- x3 y3 )
  20     swap -rot               ( x1 x2 y1 y2 )
  21     2dup * negate           ( x1 x2 y1 y2 -y1y2 )
  22     4 pick 4 pick * +       ( x1 x2 y1 y2 (x1x2-y1y2))
  23     4 roll 2 roll *         ( x2 y1 (x1x2-y1y2) x1y2 )
  24     3 roll 3 roll * +       ( (x1x2-y1y2) (x1y2+x2y1) )
  25 ;
  26
  27 : c+ ( x1 y1 x2 y2 -- x3 y3 )
  28     rot +
  29     -rot +
  30     swap
  31 ;
  32
  33 : csq 2dup c* ;
  34
  35
  36 : conj ( x y -- x -y )
  37     negate
  38 ;
  39
  40 : cmagsq ( x1 y1 -- mag )
  41     2dup conj c* +
  42 ;
  43
  44 ( --- Mandelbrot set calculations  --- )
  45
  46 : iterate ( cr ci zr zi -- cr ci z'r z'i )
  47     2over 2swap csq c+
  48 ;
  49
  50 100 value maxiter
  51
  52 : inSet? ( cr ci -- res )
  53     0 0     ( z_0 = 0 )
  54
  55     true    ( flag indicating set membership )
  56     maxiter 0 do
  57         drop
  58
  59         iterate
  60         2dup cmagsq
  61         4 0 >scaled > if
  62             false ( not in set )
  63             leave
  64         then
  65
  66         true ( maybe in set )
  67     loop
  68
  69     ( Clear z and c, leaving set membership flag )
  70     -rot 2drop -rot 2drop
  71 ;
  72
  73 100 value xsteps
  74 30 value ysteps
  75
  76 ( Draw the Mandelbrot Set!)
  77 : mandeldraw ( x1 y1 x2 y2 -- )
  78
  79     cr
  80
  81     0 pick 3 pick - ysteps /
  82     1 pick 4 pick do
  83
  84         2 pick 5 pick - xsteps /
  85         3 pick 6 pick do
  86
  87             i j inSet? if
  88                 42 emit
  89             else
  90                 space
  91             then
  92
  93         dup +loop
  94         drop
  95
  96         cr
  97
  98     dup +loop
  99     drop
 100 ;
 101
 102 ( Clean up - hide non-standard multiplication def. )
 103 hide *
 104
 105 ( Default picture )
 106 : mandel
 107         -2 0 >scaled -1 0 >scaled 0 5 >scaled 1 0 >scaled
 108         mandeldraw
 109 ;
 110
 111 CR .( Enter 'mandel' to draw the Mandelbrot Set.)