/*

 b batch file for maxima
 save as e.mac

maxima
batch("e.mac");


There are 2 complex planes :
* w-plane 
* z-plane 


z is a pont of z-planes
w is a point of w-plane

1. draw curves on z plane 
2. draw images of z-curves  under h  on w plane               
3. draw images of w-curves under hi on z plane 


http://math.stackexchange.com/questions/114733/mapping-half-plane-to-unit-disk

Adam Majewski
 
*/

kill(all);
remvalue(all);
display2d:false$







/* definitions of  functions  */


/*  map from z to w : w=h(z) */
h(z):=if z=0.0 then -1000 else rectform(float(log(z))); /* map from z plane to w plane : w = h(z) */


hi(w):= rectform(float(exp(w))); /* map from w plane to z plane : z = g(w) */


/* converts complex number into list for draw package  */
draw_format(z):= [float(realpart(z)),float(imagpart(z))];


GiveCirclePoint(t) := rectform(float(R*%e^(%i*t*2*%pi))); /* gives point of unit circle for angle t in turns */
GiveZRayPoint(R) := rectform(float(R*%e^(%i*tRay*2*%pi))); /* gives point of external ray for radius R and angle tRay in turns */ 



compile(all);


R_max: 5;
R_min: 0;
dR: R_max - R_min;
NumberOfRays: 10;
iMax: 100; /* number of points to draw */





print (" ============== compute ===============  ")$

/* curve is a list of points joned by lines */

/* f_0 plane = z-plane */
circle_angles: makelist(i/(10*iMax), i, 0, 10*iMax-1)$ /* more angles = more points */


/* Unit circle */
R: 1;
UnitCirclePointsZ: map(GiveCirclePoint, circle_angles)$
/* circle with radius 1/2 */
R: 0.5;
HalfUnitCirclePointsZ: map(GiveCirclePoint, circle_angles)$



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

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









 
 /* --------- map from z to w plane using h  -------------------------------*/

 /* w-curves = images of z-circles */
 WCirclesPoints: map (h, ZCirclesPoints)$
 UnitCirclePointsW : map(h,UnitCirclePointsZ)$
 HalfUnitCirclePointsW : map (h, HalfUnitCirclePointsZ)$
 /* w-curves = images of z-rays */
   WRaysPoints: map (h, ZRaysPoints )$

  /* --------- map from w to z plane using hi  -------------------------------*/
  
  /* z-curves = images of w-curves */
 ZiCirclesPoints: map (hi, WCirclesPoints)$
 UnitCirclePointsZi : map(hi,UnitCirclePointsW)$
 HalfUnitCirclePointsZi : map (hi, HalfUnitCirclePointsW)$
 /* w-curves = images of z-rays */
   ZiRaysPoints: map (hi, WRaysPoints )$
  
  

 

print ("--------------  convert lists of complex points to draw format lists ----------- ")$

 /* z plane */
 UnitCirclePointsZ : map (draw_format, UnitCirclePointsZ)$
 HalfUnitCirclePointsZ : map ( draw_format, HalfUnitCirclePointsZ)$
 ZCirclesPoints:map(draw_format,ZCirclesPoints)$ 
 ZRaysPoints:map(draw_format,ZRaysPoints)$ 
 
 /* w plane */
 WCirclesPoints:map(draw_format, WCirclesPoints)$ 
 WRaysPoints:map(draw_format, WRaysPoints)$ 
 UnitCirclePointsW : map (draw_format, UnitCirclePointsW)$
 HalfUnitCirclePointsW : map(draw_format, HalfUnitCirclePointsW)$

 /* z plane inverted from w */
 ZiCirclesPoints:map(draw_format, ZiCirclesPoints)$ 
 ZiRaysPoints:map(draw_format, ZiRaysPoints)$ 
 UnitCirclePointsZi : map (draw_format, UnitCirclePointsZi)$
 HalfUnitCirclePointsZi : map(draw_format, HalfUnitCirclePointsZi)$
 

print (" ------------------  draw ------------------------------------------------------ ")$

path:""$ /* pwd  */
 FileName:"j"$ /* without extension which is the terminal name */

 load(draw); /* Mario Rodríguez Riotorto   http://www.telefonica.net/web2/biomates/maxima/gpdraw/index.html */
 draw(
  terminal  = 'png,
  file_name = concat(path,FileName),
  columns  = 3,
  dimensions=[1500,500], /*  x = y*columns  */
  
  gr2d(title = " z plane ",
  /*
   yrange = [-3,3],
   xrange = [-3,3],
   */
   points_joined =true,
   grid = false,
   
   point_size    = 0.2,
   point_type = filled_circle,
   color 	= green,
   points(ZCirclesPoints),
   color         = red,
   points(UnitCirclePointsZ),
   color = yellow,
   points(HalfUnitCirclePointsZ),
   color = blue,
   points(ZRaysPoints)
   
  ),

  gr2d(
   title = " w plane: w=log(z)",
   yrange = [-2.0,2.0],
   xrange = [-2.0,2.0],
   grid = false,
   xaxis       = false,
   
   points_joined =true,

   color         = red,
   point_size    = 0.2,
   point_type = filled_circle,
   color 	= green,
   points(WCirclesPoints),
   color         = red,
   points(UnitCirclePointsW),
   color = yellow,
   points(HalfUnitCirclePointsW),
   
   color = blue,
   points(WRaysPoints)),
   
   
   
   gr2d(title = "z plane: z = exp(w) ",
  /*
   yrange = [-3,3],
   xrange = [-3,3],
   */
   points_joined =true,
   grid = false,
   
   point_size    = 0.2,
   point_type = filled_circle,
   color 	= green,
   points(ZiCirclesPoints),
   color         = red,
   points(UnitCirclePointsZi),
   color = yellow,
   points(HalfUnitCirclePointsZi),
   color = blue,
   points(ZiRaysPoints))
   
   
   
   
 );