Hi,
I have changed code. It works better ( for biger iMax gets bind stack
overflow).
Output images are slightly diffrent then in :
http://www.iec.csic.es/~miguel/Preprint4.ps
I don't know why ( now ).
Regards
Adam
================== code ====================================
/* compute c for given radius, angle and l ( lambda=multiplier of fixed
point) */
epsilon:0.001;
radius:1+epsilon;
angle:(3-sqrt(5))/2; /* in radians */
l:radius * exp(%i*angle*2*%pi);
c:l/2 - (l*l)/4;
trigsimp(%);
rectform(%);
c:float(%), numer;
/* define function ( map) for dynamical system z(n+1)=f(zn,c) */
f(z,c):=z*z+c;
/* maximal number of iterations */
iMax:200; /* to big couses bind stack overflow */
EscapeRadius:10;
/* define z-plane ( dynamical ) */
zxMin:-0.8;
zxMax:0.2;
zyMin:-0.2;
zyMax:0.8;
/* resolution is proportional to number of details and time of drawing */
iXmax:1000;
iYmax:1000;
/* compute critical point */
zcr:rhs(solve(diff(f(z,c),z,1)));
/* save critical point to 2 lists */
xcr:makelist (realpart(zcr), i, 1, 1); /* list of re(z) */
ycr:makelist (imagpart(zcr), i, 1, 1); /* list of im(z) */
/* ------------------- compute forward orbit of critical point ----------*/
z:zcr; /* first point */
orbit:[z];
for i:1 thru iMax step 1 do
block
(
z:f(z,c),
if abs(z)>EscapeRadius then return(i) else orbit:endcons(z,orbit)
);
/*-------------- save orbit to draw it later on the screen
----------------------------- */
/* save the z values to 2 lists */
xx:makelist (realpart(zcr), i, 1, 1); /* list of re(z) */
yy:makelist (imagpart(zcr), i, 1, 1); /* list of im(z) */
for i:2 thru length(orbit) step 1 do
block
(
xx:cons(realpart(orbit[i]),xx),
yy:cons(imagpart(orbit[i]),yy)
);
/* draw reversed orbit of beta using draw package */
load(draw);
draw2d(
terminal = 'screen,
pic_width = iXmax,
pic_height = iYmax,
yrange = [zyMin,zyMax],
xrange = [zxMin,zyMax],
title= concat(" dynamical plane for f(z,c):=z*z+",string(c)),
xlabel = "Z.re ",
ylabel = "Z.im",
point_type = filled_circle,
points_joined = false,
color =red,
point_size = 0.5,
key = "critical orbit",
points(xx,yy),
key = "critical point",
color =blue,
points_joined = false,
points(xcr,ycr)
);