kill(all) $
f(z,c):=z*z+c $
critical_orbit(c,iMax):=
block(
/* iMax : Number of iterations , for most cases iMax:100 */
ER:2.0,
z:0+0*%i, /* first point = critical point */
orbit:[z],
/* compute forward orbit */
i:0,
if (abs(z)>ER) then return(orbit),
loop,
z:rectform(f(z,c)),
orbit:endcons(z,orbit),
i:i+1,
if ((abs(z)<ER) and (i<iMax)) then go(loop),
return(orbit)
)$
/* find fixed point alfa */
GiveFixed(c):= float(rectform((1-sqrt(1-4*c))/2))$
GiveDistanceFromCenter(z):= abs(z-zf)$
GiveInnerRadiusOfOrbit(orbit):=lmin(map(GiveDistanceFromCenter,orbit))$
GiveOuterRadiusOfOrbit(orbit):=lmax(map(GiveDistanceFromCenter,orbit))$
/* =========== main =========*/
c:-.5867879078859505*%i-.3905408691260131; /* Golden Mean*/
zf:GiveFixed(c);
zfx:realpart(zf)$
zfy:imagpart(zf)$
iXMax:200;
NrPoints:10*iXMax;
orbit:critical_orbit(c,NrPoints)$
innerRadius: GiveInnerRadiusOfOrbit(orbit) ;
ir2 : innerRadius * innerRadius $
outerRadius: GiveOuterRadiusOfOrbit(orbit) ;
or2 : outerRadius * outerRadius $
load(draw);
draw2d(
title= concat("Critical orbit, inner and outer circle for fc(z)=z*z", string(c)),
user_preamble = "set size ratio 1; set key outside right; set key box ", /* */
file_name = "croGM2d",
terminal = png,
yrange = [-0.8,0.2],
xrange = [-0.8,0.2],
dimensions = [800,800],
xlabel = "Z.re ",
ylabel = "Z.im",
point_type = filled_circle,
points_joined = false,
color =blue,
point_size = 0.9,
key = "center ",
points([[realpart(zf),imagpart(zf)]]),
key = "critical orbit",
color =red,
point_size = 0.5,
points(map(realpart,orbit),map(imagpart,orbit)),
ip_grid = [400,400], /* Number of initial grid points in implicit plots */
key = "inner circle",
color = green,
implicit( (x-zfx)^2+(y-zfy)^2 = ir2 , x,-2,2,y,-2,2),
key = "outer circle",
color = black,
implicit( (x-zfx)^2+(y-zfy)^2 = or2 , x,-2,2,y,-2,2),
color = magenta ,
point_size = 0.9,
key = "crital point ",
points([[0,0]])
);