Hi,
Here is the program which IMHO needs more memory
Try change jMax to 50 or 100 to see effect.
( result image : http://commons.wikimedia.org/wiki/File:Jung50e.png
theory : http://fraktal.republika.pl/mset_jungreis.html
or
http://en.wikipedia.org/wiki/External_ray
Regards
Adam
-====================== code ===================
/* batch file for maxima
uses :
- symmetry around horizontal ( 0X ) axis
- Psi_M function to map conjugate plane to parameter plane
- jungreis algorithm to
time for jMax:30; iMax:200; iMaxBig:400; is 1191 sec
*/
start:elapsed_run_time ();
jMax:30; /* precision = proportional to details and time of
computations */
iMax:200; /* number of points to draw */
iMaxBig:400;
/* computes b coefficient of Jungreis function*/
betaF[n,m]:=block
(
[nnn:2^(n+1)-1],
if m=0
then 1.0
else if ((n>0) and (m < nnn))
then 0.0
else (betaF[n+1,m]-
sum(betaF[n,k]*betaF[n,m-k],k,nnn,m-nnn)-betaF[0,m-nnn])/2.0
)$
b[m]:=betaF[0,m+1]$
/* -------------------------------*/
/* Power of w to j */
wn[w,j]:= if j=0 then 1 else w*wn[w,j-1]$
/* ---------Jungreis function ; c = Psi_M(w)
----------------------------- */
Psi_M(w):=w + sum(b[j]/wn[w,j],j,0,jMax)$
/* exponential for of complex number with angle in turns */
GiveCirclePoint(t):=R*%e^(%i*t*2*%pi)$ /* gives point of unit circle for
angle t in turns */
GiveWRayPoint(R):=R*%e^(%i*tRay*2*%pi)$ /* gives point of external ray
for radius R and angle tRay in turns */
NmbrOfCurves:9;
/* coordinate of w-plane not c-plane */
R_max:1.5;
R_min:1;
dR:(R_max-R_min)/NmbrOfCurves; /* for circles */
dRR:(R_max-R_min)/iMax; /* for rays */
/* --------------------------------------f_0 plane = w-plane
-----------------------------------------*/
/*-------------- unit circle ------------*/
R:1;
circle_angles:makelist(i/iMax,i,0,iMax/2)$
CirclePoints:map(GiveCirclePoint,circle_angles)$
/* ---------------external circles ---------*/
circle_radii: makelist(R_min+i*dR,i,1,NmbrOfCurves)$
circle_angles:makelist(i/iMaxBig,i,0,iMaxBig/2)$
WCirclesPoints:[]$
for R in circle_radii do
WCirclesPoints:append(WCirclesPoints,map(GiveCirclePoint,circle_angles))$
/* -------------- external w rays -------------*/
ray_radii:makelist(R_min+dRR*i,i,1,iMax);
ray_angles:[0,1/3,1/7 , 2/7 ,3/7 ]; /* list of angles < 1/2 of root
points */
WRaysPoints:[];
for tRay in ray_angles do
WRaysPoints:append(WRaysPoints,map(GiveWRayPoint,ray_radii));
/* -------------------------parameter plane = c plane
-----------------------------------*/
MPoints:map(Psi_M,CirclePoints); /* Mandelbrot set points */
CRaysPoints:map(Psi_M,WRaysPoints); /* external z rays */
Equipotentials:map(Psi_M,WCirclesPoints);
/* add points below horizontal axis */
for w in CirclePoints do CirclePoints:cons(conjugate(w),CirclePoints);
for w in WRaysPoints do WRaysPoints:cons(conjugate(w),WRaysPoints);
for w in WCirclesPoints do WCirclesPoints:cons(conjugate(w),WCirclesPoints);
for c in MPoints do MPoints:cons(conjugate(c),MPoints);
for c in CRaysPoints do CRaysPoints:cons(conjugate(c),CRaysPoints);
for c in Equipotentials do Equipotentials:cons(conjugate(c),Equipotentials);
/* time */
stop:elapsed_run_time ();
time:fix(stop-start);
/* ---------------- draw
*--------------------------------------------------------------------------*/
load(draw); /* Mario Rodr?guez Riotorto
http://www.telefonica.net/web2/biomates/maxima/gpdraw/index.html */
draw(file_name = "jung30es_2",
pic_width=1000,
pic_height= 500,
terminal = 'png,
columns = 2,
gr2d(title = " unit circle with external rays & circles ",
point_type = filled_circle,
points_joined =true,
point_size = 0.34,
color = red,
points(map(realpart, CirclePoints),map(imagpart, CirclePoints)),
points_joined =false,
color = green,
points(map(realpart, WCirclesPoints),map(imagpart, WCirclesPoints)),
color = black,
points(map(realpart, WRaysPoints),map(imagpart, WRaysPoints))
),
gr2d(title = "Parameter plane : Image under Psi_M(w) ",
points_joined =true,
point_type = filled_circle,
point_size =0.34,
color = blue,
points(map(realpart, MPoints),map(imagpart, MPoints)),
points_joined =false,
color = green,
points(map(realpart, Equipotentials),map(imagpart, Equipotentials)),
color = black,
points(map(realpart, CRaysPoints),map(imagpart, CRaysPoints))
)
);