Improved error reporting.
[forth.jl.git] / examples / mandelbrot.4th
index c75ebd0..9759a39 100644 (file)
 ( --- Complex arithmetic --- )
 
 ( Location of floating point. )
-: precision 10000 ;
+10000 value precision
 
-: >scaled precision 10 / * swap precision * + ;
+: sign dup abs / ;
+
+: >scaled
+    precision 10 / * over
+    ?dup 0<> if 
+        sign *
+    then
+    swap precision * +
+;
 
 ( Redefine multiplication.  Yay forth! )
 : * precision */ ;
 
 : c* ( x1 y1 x2 y2 -- x3 y3 )
-        swap -rot               ( x1 x2 y1 y2 )
-        2dup * negate           ( x1 x2 y1 y2 -y1y2 )
-        4 pick 4 pick * +       ( x1 x2 y1 y2 (x1x2-y1y2))
-        4 roll 2 roll *         ( x2 y1 (x1x2-y1y2) x1y2 )
-        3 roll 3 roll * +       ( (x1x2-y1y2) (x1y2+x2y1) )
+    swap -rot               ( x1 x2 y1 y2 )
+    2dup * negate           ( x1 x2 y1 y2 -y1y2 )
+    4 pick 4 pick * +       ( x1 x2 y1 y2 (x1x2-y1y2))
+    4 roll 2 roll *         ( x2 y1 (x1x2-y1y2) x1y2 )
+    3 roll 3 roll * +       ( (x1x2-y1y2) (x1y2+x2y1) )
 ;
 
 : c+ ( x1 y1 x2 y2 -- x3 y3 )
-        rot +
-        -rot +
-        swap
+    rot +
+    -rot +
+    swap
 ;
 
 : csq 2dup c* ;
 
+
+: conj ( x y -- x -y )
+    negate
+;
+
 : cmagsq ( x1 y1 -- mag )
-        csq abs
+    2dup conj c* +
 ;
 
 ( --- Mandelbrot set calculations  --- )
 
 : iterate ( cr ci zr zi -- cr ci z'r z'i )
-        csq c+
+    2over 2swap csq c+
 ;
 
+100 value maxiter
+
 : inSet? ( cr ci -- res )
+    0 0     ( z_0 = 0 )
     
-    100 0 do
+    true    ( flag indicating set membership )
+    maxiter 0 do
+        drop
 
-        2swap 2dup 5 roll 5 roll
         iterate
         2dup cmagsq
-        100 > if
+        4 0 >scaled > if
+            false ( not in set )
             leave
         then
 
+        true ( maybe in set )
     loop
+
+    ( Clear z and c, leaving set membership flag ) 
+    -rot 2drop -rot 2drop
 ;
 
+100 value xsteps
+30 value ysteps
+
 ( Draw the Mandelbrot Set!)
-: mandel ( x1 y1 x2 y2 -- )
+: mandeldraw ( x1 y1 x2 y2 -- )
+
+    cr
 
+    0 pick 3 pick - ysteps /
+    1 pick 4 pick do
+
+        2 pick 5 pick - xsteps /
+        3 pick 6 pick do
+
+            i j inSet? if
+                42 emit
+            else
+                space
+            then
+
+        dup +loop
+        drop
+
+        cr
+
+    dup +loop
+    drop
 ;
 
 ( Clean up - hide non-standard multiplication def. )
 hide *
+
+( Default picture )
+: mandel
+        -2 0 >scaled -1 0 >scaled 0 5 >scaled 1 0 >scaled
+        mandeldraw
+;
+
+CR .( Enter 'mandel' to draw the Mandelbrot Set.)