R_max: 5;
R_min: 1;
dR: R_max - R_min;
psi(w) := w+1/w;
NmbrOfRays: 10;
iMax: 100; /* number of points to draw */
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 */ 

/* f_0 plane = W-plane */
/* Unit circle */
R: 1;
circle_angles: makelist(i/(10*iMax), i, 0, 10*iMax-1); /* more angles = more points */
CirclePoints: map(GiveCirclePoint, circle_angles);

/* External circles */
circle_radii: makelist(R_min+i, i, 1, dR);
WCirclesPoints: [];
for R in circle_radii do 
	WCirclesPoints: append(WCirclesPoints, map(GiveCirclePoint, circle_angles));

/* External W rays */
ray_radii: makelist(R_min+dR*i/iMax, i, 0, iMax);
ray_angles: makelist(i/NmbrOfRays, i, 0, NmbrOfRays-1);
WRaysPoints: [];
for tRay in ray_angles do 
	WRaysPoints: append(WRaysPoints, map(GiveWRayPoint, ray_radii));


/* f_c plane = Z plane = dynamic plane */
/* external Z rays */
ZRaysPoints: map(psi, WRaysPoints);

/* Julia set points */
JuliaPoints: map(psi, CirclePoints);
Equipotentials: map(psi, WCirclesPoints);


/* Mario Rodríguez Riotorto (http://www.telefonica.net/web2/biomates/maxima/gpdraw/index.html) */
load(draw);
draw(
	file_name = "erays",
	pic_width = 1000, 
	pic_height = 500,
	terminal = 'svg,
	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 = black,
		points(map(realpart, WRaysPoints), map(imagpart, WRaysPoints)),
		points(map(realpart, WCirclesPoints), map(imagpart, WCirclesPoints))
	),
	gr2d(
		title = "Image under psi(w):=w+1/w; ",
		points_joined = true,
		point_type = filled_circle,
		point_size = 0.34,
		color = blue,
		points(map(realpart, JuliaPoints),map(imagpart, JuliaPoints)),
		points_joined = false,
		color = black,
		points(map(realpart, ZRaysPoints),map(imagpart, ZRaysPoints)),
		points(map(realpart, Equipotentials),map(imagpart, Equipotentials))
	) 
);