> # Program  to make an animated plot of Julia sets

>      Based on a program from "Dynamical Systems With Applications Using MAPLE" , by 
      Stephen Lynch, Birkhauser 2001, Chapter 15. This algorithm implements the Chaos Game.
> ##################################
      This program produces an animated plot of a sequence of Julia sets, parameterized by c.
      Store the plot as an animated gif file and you can insert it inot a webpage.
      Time taken for calculating: M=25, k=12;  about 30 minutes ( Celeron 500) .

      Here a Julia set will be stored in arrays holding the x -coordinates and y-coordinates of the 
      points (x1 and y1). 
> 
> 
> M:=38:     #    number of interpolating steps between initial c and final c (= # of steps in movie). Default = 38
> k:=12:iter:=2^k:  # iter is the number of points in the Julia set. Defualt = 12
> ao:=0.1:bo:=0.8:     #  initial c=(ao, bo)
> af:=-0.3:bf:=0.8:    #  final c=(af, bf)
> da:=(af-ao)/M:db:=(bf-bo)/M:   # da and db are the step sizes of the increments (the smaller the 
                                                                               #  they are the smoother the movie appears; 0.06 looks fine to me)
> x1:=array(0..M,0..iter+1):y1:=array(0..M,0..iter+1):  # these arrays hold the Julia sets,
                                # first index is the value of c and the second index are the x-coordinates (y-coordinate) of the points.
> 
> for j from 0 to M do          # beginning of the outer for loop (for each j compute a Julia set)
> a:=ao+j*da:b:=bo+j*db:      #  increment a by da and b by db at each step of the outer loop
> die:=rand(0..1):  #generates random numbers 0 or 1
> # Determine the unstable fixed point 
> x1[j,0]:=Re(0.5+sqrt(0.25-(a+I*b))):
> y1[j,0]:=Im(0.5+sqrt(0.25-(a+I*b))):
> 2*abs(x1[j,0]+I*y1[j,0]);
> for i from 0 to iter do     #   beginning of the inner for loop to draw the points of the Julia set
> x:=x1[j,i]:y:=y1[j,i]:
> u:=sqrt((x-a)^2+(y-b)^2)/2:v:=(x-a)/2:
> u1:=evalf(sqrt(u+v)):v1:=evalf(sqrt(u-v)):
> x1[j,i+1]:=u1:y1[j,i+1]:=v1:if y1[j,i]<b then y1[j,i+1]:=-y1[j,i+1]:fi:
> die();
> if (die()=0) then x1[j,i+1]:=-u1:y1[j,i+1]:=-y1[j,i+1]:fi:od:
              #    end of inner for loop
> od:   #  end of outer for loop
  Plotting
> with(plots,animate):animate([x1[t,m],y1[t,m],m=0..iter],t=0..M,numpoints=iter+1,color= black,style=point,symbol=point,frames=M+1);
     This plots the array containing the Julia sets (parameterized by the first index, which is the value of c). It is treated as a parameterized 
     family of curves, the parameter t running from 0 (initial c) to M (final c). Once this plot is made you can save it using the export option in the 
    plotting menu; save as a gif file. 
>