############################################################################
#
#	File:     julia1.icn
#
#	Subject:  Program to display the Julia set
#
#	Author:   Ralph E. Griswold
#
#	Date:     June 17, 1994
#
############################################################################
#
#   This file is in the public domain.
#
############################################################################
#
#  This is a barebones version of a display of the Julia set.  It
#  has deliberately been left simple and free of options so that the
#  basic idea is clear and so that it can be used as the basis of
#  more capable versions.
#
#  This program is based on material given in "Chaos, Fractals,
#  and Dynamics", Robert L. Devaney, Addison-Wesley, 1990.
#
#  The point in the complex plane for which the Julia set is computed
#  is given on the command line, as in
#
#	julia1 .360284 .100376
#
#  which displays the Julia set for the complex number .360284 + .100376i.
#
############################################################################
#
#  Requires:  Version 9 graphics
#
############################################################################
#
#  Links:  wopen
#
############################################################################

link wopen

procedure main(argl)
   local c1, c2, extent, half, quarter, m, n, x0, y0, x, y
   local x1, y1, i, z

   c1 := real(argl[1]) | -1.0		# default is -1.0 + 0.0i
   c2 := real(argl[2]) | 0.0

   extent := 200
   half := 200 / 2
   quarter := real(extent) / 4

   WOpen("label=julia", "height=" || extent, "width=" || extent) |
      stop("*** cannot open window")

   every m := 0 to extent do {
      x0 := -2 + m / quarter
      every n := 0 to half do {
         y0 := 2 - n / quarter
         x := x0
         y := y0
         every i := 1 to 20 do {	# compute orbit
            x1 := x ^ 2 - y ^ 2 + c1
            y1 := 2 * x * y + c2
            x := x1
            y := y1
            z := x ^ 2 + y ^ 2
            if z > 4 then break next	# if escaping, forget it
            }
         DrawPoint(m, n)
         DrawPoint(extent - m, extent - n)
         }
      }

   Event()

end