# File: fraktale.txt
# All procedures were written for Maple V Release 3.  
# 
# Original file location: 
# http://www.math.utsa.edu/mirrors/maple/maplev/fractals/fraktale.ms
# 
# ------------------------------------------------------------------------------------------------
# 
# Note: Make sure that Digits <= 10 in order to speed up computations.
# 
# Credits: 
# John Oprea, oprea@math.csuohio.edu, implemetented the J procedure to pass initial
# points and the anonymous functions to the fractal procedures. His procedures also
# caused me to add further fractals.
# 
# The escape-time fractal procedures are based on a Maple V algorithm taken from the book 
# 'Maple V - Programming Guide' by M.B. Monagan, K.O. Geddes, G. Labahn and S. Vorkoetter, 
# Springer Verlag.
# 
# Remark on the first procedure 'mandelbrot' by John Oprea:
# The modification on mandelbrot uses the fact that Maple can color things according to a function. 
# Also, plotting the zero function keeps the graph in a plane.  Increasing the values of 'grid' in the 
# plot3d command results in  better resoultion of the fractal. However, this means longer  
# computation time.  'style=patchnogrid' draws the surface with shaded rectangular patches  with 
# no wireframe grid. This is the best 'style' option for this routine.  The procedure mandelbrot returns 
# the value m. The higher m the more  likely the probability that the point (x, y) is part of the 
# Mandelbrot  set. Certainly m depends of the number of iterations (here: 30). It  is possible that a 
# point satisfies the condition abs(z)<2 after n  iterations but is outside the range after n+1 iterations. 
# So to make sure that (x, y) is part of the set you would have to do infinite iterations.  m is used in the 
# 'color=mandelbrot' option to visualize the set:  the red area - the Mandelbrot lake - indicates that
# it is part of the set.  
> restart:


> mandelbrot:=proc(x, y)
>    local z, c, m;
>    c:=evalf(x+y*I);
>    z:=evalf(x+y*I);
>    for m from 0 to 30 while abs(z)<2 do
>       z:=z^2+c
>    od;
>    m
> end:


> plot3d(0, -2 .. 0.7, -1.2 .. 1.2, orientation=[-90,0], grid=[250, 250], style=patchnogrid, 
> scaling=constrained, color=mandelbrot);


> plot3d(0, -1.3 .. -1.25, 0.35 .. 0.4, orientation=[-90,0], grid=[100, 100], style=patchnogrid, scaling=con
> strained, color=mandelbrot);
> restart: 


> mandelbrot:=proc(x, y)
>    local z, m;
>    z:=evalf(x+y*I);
>    for m from 0 to 25 while abs(z)<2 do
>       z:=z^2+(x+y*I)
>    od;
>    m
> end:


> plot3d(mandelbrot, -2 .. 0.5, -1.2 .. 1.2, grid=[50, 50]);


> restart:


# Julia sets modified by John Oprea, oprea@math.csuohio.edu  
> julia := proc(c, x, y)
>    local z, m;
>    z := evalf(x+y*I);
>    for m from 0 to 30 while abs(z) < 3 do
>       z := z^2+c
>    od;
>    m
> end:


> J := proc(d)
>    global phonyvar;
>    phonyvar := d;
>    (x, y) -> julia(phonyvar, x, y)
> end:


> plot3d(0, -2 .. 2, -1.3 ..1.3, style=patchnogrid, orientation=[-90,0], grid=[150, 150],
> scaling=constrained, color=J(-1.25));


> plot3d(0, -1.9 .. 1.9, -1.9 ..1.9, style=patchnogrid, orientation=[-90,0], grid=[100, 100],
> scaling=constrained, color=J(-.74543+.11301*I));
> plot3d(0, -1.9 .. 1.9, -1.2 ..1.2, style=patchnogrid, orientation=[-90,0], grid=[150, 150],
> scaling=constrained, color=J(-.9+.12*I));
> restart:


> # LAMBDAFN
> lambdafn_sin:=proc(x, y)
>    local z, m, p1;
>    z:=evalf(x+y*I);
>    p1:=evalf(1+I*0.4);
>    for m from 0 to 100 while abs(z)<4 do
>       z:=sin(z)*p1
>    od;
>    m
> end:


> plot3d(0, -4 .. 4, -3 .. 3, orientation=[-90,0], grid=[50, 50], style=patchnogrid, 
> scaling=constrained, color=lambdafn_sin);


> plot3d(0, -1.5 .. 1.5, -1.1 .. 1.1, orientation=[-90,0], grid=[250, 250], style=patchnogrid, 
> scaling=constrained, color=lambdafn_sin);
> restart: 


> # BIOMORPH2
> # May 22, 1996
> biomorph2:=proc(x, y)
>    local z, m, c;
>    z:=evalf(x+y*I);
>    c:=z;
>    for m from 0 to 100 while abs(z)<4 do
>       z:=cos(z)*z^2+c
>    od;
>    m
> end:


> plot3d(0, -2 .. 2, -1.5 .. 1.5, orientation=[-90,0], grid=[20, 20], style=patchnogrid, 
> scaling=constrained, color=biomorph2);


> restart: 
> # NEWTON 
> # Visualization of Newton's method to calculate the complex solution of z^n-1=0 by 
> # the iteration z(j+1) = z(j) - fn(z)/diff(fn(z)) with fn(z) := z^n-1, and n > 2 the polynomial degree
> # of fn (also real values for n are allowed).
> solve(z^3-1=0, z);
> # The graph shows the starting points from which the method does or does not 
> # converge to one of the above roots. The documentation of FRACTINT explains the result: 
> # "The coloring of the plot shows the "basins of attraction" for each root of the polynomial 
> # -- i.e., an initial guess within any area of a given color would lead you to one of the roots. 
> # As you can see, things get [fractal] along certain radial lines or "spokes," those being the lines
> # between actual roots." 
> #
> # The fractal properties of Netwon's method to calculate complex roots were first
> # discovered by John Hubbard of Orsay, France.
> newton:=proc(x, y)
>    local z, m;
>    z:=evalf(x+y*I);
>    for m from 0 to 50 while abs(z^3-1) >= 0.001 do
>       z:=z-(z^3-1)/(3*z^2)
>    od;
>    m
> end:
> plot3d(0, -2 .. 2, -1.5 .. 1.5, orientation=[-90,0], grid=[250, 250], style=patchnogrid, 
> scaling=constrained, color=newton);
> restart:
> cmplxmarkmand:=proc(x, y)
> # May 25, 1996
> # The original formula 'z:=z^2*c^(p-1)+c' was developed by Mark Peterson.
>    local c, m, z;
>    c:=evalf(x+y*I);
>    z:=evalf(x+y*I);
>    for m from 0 to 30 while abs(z)<2 do
>       z:=z^2*c^0.1+c
>    od;
>    ln(m)
> end:


> plot3d(0, -2 .. 0.7, -1.2 .. 1.2, orientation=[-90,0], grid=[250, 250], style=patchnogrid, 
> scaling=constrained, color=cmplxmarkmand);
> restart:


> lambda:=proc(x, y)
> # Lambda
> # May 26, 1996
> # 
> # Complex version of the population growth equation x(j+1)=r*x(j)*(1-x(j))
> # developed by P. F. Verhulst in 1845
> # Lambda is part of the Julia family
>    local c, z, m;
>    c:=0.85+I*0.6;
>    z:=evalf(x+y*I);
>    for m from 0 to 30 while abs(z)<3 do
>       z:=c*z*(1-z)
>    od;
>    m
> end:


> plot3d(0, -1.5 .. 2.5, -1.5 ..1.5, style=patchnogrid, orientation=[-90,0], grid=[25, 25], 
> scaling=constrained, color=lambda);


> restart:

> ifs:=proc(imax)
> # Iterated Functions Systems
> # original program 'Chaos-Spiel fuer ein Farnblatt' written by Heinz-Otto Peitgen,
> # Hartmut Juergens and Dietmar Saupe in BASIC,
> # published in: 'Bausteine des Chaos - Fraktale', Springer Verlag/Klett-Cotta, 
> # p. 415, written by the authors mentioned above
> # June 06, 1997
>    local e1, e2, e3, e4, f1, f2, f3, f4, x, y, xn, yn, z, pts;
>    pts:=NULL;
>    e1:=0.5; e2:=0.57; e3:=0.408; e4:=0.1075;
>    f1:=0; f2:=-0.036; f3:=0.0893; f4:=0.27;
>    x:=e1;
>    y:=0;
>    to imax do
>       z:=rand()/1e12;
>       if z<=0.02 then 
>          xn:=e1;
>          yn:=0.27*y+f1
>       elif z<=0.17 then 
>          xn:=-0.139*x+0.263*y+e2;
>          yn:=0.246*x+0.224*y+f2
>       elif z<=0.3 then 
>          xn:=0.17*x-0.215*y+e3;
>          yn:=0.222*x+0.176*y+f3
>       else
>          xn:=0.781*x+0.034*y+e4;
>          yn:=-0.032*x+0.739*y+f4
>       fi;
>       pts:=pts, [xn, yn];
>       x:=xn;
>       y:=yn
>    od;
>    [pts]
> end:


> plot(ifs(12500), axes=NONE, style=point, symbol=POINT, scaling=constrained);


> restart:


> julfn_zsqrd:=proc(x, y)
> # Based on a Maple V algorithm taken from the book 
> # 'Maple V - Programming Guide' by M.B. Monagan, 
> # K.O. Geddes, G. Labahn and S. Vorkoetter, 
> # Springer Verlag, modified by John Oprea
> # modified by Alexander F. Walz
> # May 25, 1996
>    local c, z, m;
>    c:=evalf(-0.5+0.5*I);
>    z:=evalf(x+y*I);
>    for m from 0 to 30 while abs(z)<3 do
>       z:=sin(z)+z^2+c
>    od;
>    m
> end:


> plot3d(0, -1.91 .. 1.37, -1.24 ..1.21, style=patchnogrid, orientation=[-90,0], grid=[250, 250], scaling=co
> nstrained, color=julfn_zsqrd);
> restart: 


> mandelbrot_fast:=proc(x, y)
> # Based on a Maple V algorithm taken from the book 
> # 'Maple V - Programming Guide' by M.B. Monagan, 
> # K.O. Geddes, G. Labahn and S. Vorkoetter, 
> # Springer Verlag, modified by John Oprea
> # modified by Alexander F. Walz
> # May 25, 1996
>    local xn, xnold, yn, m;
>    xn:=x;
>    yn:=y;
>    for m from 0 to 100 while sqrt(evalhf(xn^2+yn^2)) < 2 do
>       xnold:=xn;   
>       xn:=evalhf(xn^2-yn^2+x);
>       yn:=evalhf(2*xnold*yn+y)
>    od;
>    m
> end:


> plot3d(0, -0.845 .. -0.794, 0.172 .. 0.21, orientation=[-90,0], grid=[250, 250], style=patchnogrid, 
>    scaling=constrained, color=mandelbrot_fast);


> restart:

> chaos:=proc(x1, y1, x2, y2, x3, y3, maxiter, seed1, seed2)
> # original procedure written by Tom Williams and modified by John Oprea
>    local randi, x, y, sx, sy, ir, i, j;
>    randi:=rand(1 .. 3):
>    x:=[x1, x2, x3]:
>    y:=[y1, y2, y3]:
>    sx:=seed1:
>    sy:=seed2:
>    ir:=[sx, sy]:
>    for i from 1 to maxiter do   
>       j:=randi():   
>       sx:=(sx+x[j])/2.:
>       sy:=(sy+y[j])/2.:
>       ir:=[op(ir), sx, sy]  
>    od:
>    plot(ir, style=point, symbol=POINT, scaling=constrained)
> end:


> chaos(0, 0, .5, 1, 1, 0, 9000, .5, .5);


> chaos(0, 0, -4, 7, 6, 2, 3000, 1, .5);
> restart:

> alexerror:=proc(x, iter)
> # written by Alexander F. Walz
> # May 29, 1996
> # March 02, 1997
>    local n, pts, xn, xnold, yn, y;
>    pts:=[0, 0];  # in case pts is not assigned any value after completion of the for n loop
>    for y from 0 to 1.1 by 0.01 do
>       xn := x;
>       yn := y;
>       for n to iter while evalhf(xn^2+yn^2) < 4 do
>          xnold:=xn;
>          xn := evalhf(xn^2-yn^2+x); 
>          yn := evalhf(2*xnold*yn+y)
>       od;
>       if n >= iter then pts := pts, [xn, yn], [xn, -yn] fi
>       # as a variant try: "if n >= iter then pts := pts, [x, y], [-x, -y] fi;" to compute the Mandelbrot set  
>    od; 
>    pts 
> end:


> i:='i': PLOT(seq(POINTS(alexerror(i/100, 500)), i=-200 .. 70), SYMBOL(POINT), axes=none);


> restart: 

> orbits:=proc(x, y, iter)
> # procedure to plot the trajectory a point traverses by being iterated 
> # through the formula 'z(j+1)=z(j)^2+c' (also called 'orbits') for the
> # Mandelbrot set
> # written by Alexander F. Walz
> # May 30, 1996
>    local a, b, xn, xnold, yn, pts;
>    pts:=NULL;
>    a:=x; 
>    b:=y; 
>    xn:=x;
>    yn:=y;
>    to iter do
>       xnold:=xn;   
>       xn:=evalf(xn^2-yn^2+a);
>       yn:=evalf(2*xnold*yn+b);
>       pts:=pts, [xn, yn]
>    od; 
>    [pts]
> end:


> # orbits1.gif   'Spiral'
> plot(orbits(-0.61143695190549, -0.36571035906672, 500), style=point, symbol=POINT, 
> axes=framed, scaling=constrained);


> # orbits2.gif 
> plot(orbits(-0.68181818351150, -0.31095156446099, 300), style=point, symbol=POINT, 
> axes=framed, scaling=constrained);
> # orbits3.gif   'Vortex'
> plot(orbits(-0.34555229917169, -0.60430224984884, 300), style=point, symbol=POINT, 
> axes=framed, scaling=constrained);
> # orbits4.gif
> plot(orbits(+0.24877810105681, +0.58083451911807, 115), style=line, 
> axes=framed, scaling=constrained);
> # orbits5.gif   'Radiation'
> plot(orbits(-0.33382209390402, -0.60821359232068, 500), style=point, symbol=POINT, 
> axes=framed, scaling=constrained);
> # orbits6.gif
> plot(orbits(-0.25953079387546, -0.62385896220803, 500), style=point, symbol=POINT, 
> axes=framed, scaling=constrained);
> # orbits7a.gif   'Curved'
> plot(orbits(-0.14222874119878, -0.64732701703906, 500), style=point, symbol=POINT, 
> axes=framed, scaling=constrained);
> # orbits7b.gif   'Coil'
> plot(orbits(-0.14222874119878, -0.64732701703906, 300), style=line, 
> axes=framed, scaling=constrained);
> # orbits8a.gif   'Spiral Galaxies'
> plot(orbits(-0.50977517291904, -0.60039090737700, 500), style=point, symbol=POINT, 
> axes=framed, scaling=constrained);
> # orbits8b.gif   'US Star'
> plot(orbits(-0.50977517291904, -0.60039090737700, 300), style=line, 
> axes=framed, scaling=constrained);
> restart:

> # LAMBDAFN_FAST
> # Using hardware coprocessor
> # Based on a Maple V algorithm taken from the book 
> # 'Maple V - Programming Guide' by M.B. Monagan, 
> # K.O. Geddes, G. Labahn and S. Vorkoetter, 
> # Springer Verlag, modified by John Oprea
> # modified by Alexander F. Walz
> # May 31, 1996
> #
> lambdafn_sinfast:=proc(x, y)
>    local xn, xnold, yn, m;
>    xn:=x;
>    yn:=y;
>    for m from 0 to 100 while evalhf(sqrt(xn^2+yn^2)) < 4 do
>       xnold:=xn;
>       xn:=evalhf(sin(xn)*cosh(yn) - cos(xn)*sinh(yn)*0.4);
>       yn:=evalhf(cos(xnold)*sinh(yn) + sin(xnold)*cosh(yn)*0.4)
>    od;
>    m
> end:


> plot3d(0, -4 .. 4, -3 .. 3, orientation=[-90,0], grid=[250, 250], style=patchnogrid, scaling=constrained, 
> color=lambdafn_sinfast);
> plot3d(0, -1.5 .. 1.5, -1.1 .. 1.1, orientation=[-90,0], grid=[25, 25], style=patchnogrid, scaling=constrai
> ned, color=lambdafn_sinfast);
> restart:

> kamtorus:=proc(angle, step, ende, points)
> # written by Alexander F. Walz
> # May 31, 1996
>    local x, xold, y, orbit, pts;
>    pts:=NULL;
>    for orbit from 0 by step to ende do
>       x:=orbit/3;
>       y:=x;
>       to points 
>         while abs(x) < 10^38 or abs(y) < 10^38 do 
>            xold:=x;
>            x := evalhf(x*cos(angle) + (x^2-y)*sin(angle));
>            y := evalhf(xold*sin(angle) - (xold^2-y)*cos(angle));
>            pts:=pts, [x, y]
>         od
>    od; 
>    [pts]
> end:


> plot(kamtorus(1.3, 0.05, 1.5, 200), style=point, symbol=POINT, axes=none, scaling=constrained);


> restart:

> koch := proc(p: numeric)
> # SNOWFLAKE
> # developed by Swedish mathematician Helge von Koch in 1904
> # algorithm taken from 'Algorithmen fuer Chaos und Fraktale'
> # by Dietmar Herrmann, Addison-Wesley
> # original program written in BASIC
> # p is the number of recursions, should be no more than value 5
>    local m, n, k, l, s, h, x, y, pts, t;
>    h := 3^(-p);
>    pts := NULL;
>    x := 0; y := 0;
>    for n from 0 to (4^p-1) do
>       m := n;
>       s := 0;
>       for l from 0 to p-1 do
>          t := irem(m, 4);
>          m := iquo(m, 4);
>          s := s+irem((t+1), 3) -1
>       od;
>       x := x+cos(Pi*s/3)*h;
>       y := y+sin(Pi*s/3)*h;
>       pts := pts, [x, y]
>    od; 
>    [pts]
> end:


> plot(koch(5), scaling=constrained, style=LINE, axes=NONE);