/*  
 draws external dynamic rays 
 R(t) = {z:arg_e(z)=t}
 using 
 z= Psi_n(w) = fc^{-n}(w^2^n) 
 there are 2 dynamic planes : 
 - f0 plane where are w points; f0(w):=w*w
 - fc plane where are z points; fc(z):=z*z+c
 */ 


kill(all);
remvalue(all);


 /* --------------------------definitions of functions ------------------------------*/
 f(z,c):=z*z+c;
 finverseplus(z,c):=sqrt(z-c);
 finverseminus(z,c):=-sqrt(z-c); 

 /*
 Square root of complex number : csqrt(x + y * i) = sqrt((r + x) / 2) + i * y / sqrt(2 * (r + x))
 gives principal value of square root : -Pi <arg<Pi
 */
 csqrt(z):=
 block(
  [t,re,im],
  t:abs(z)+realpart(z),
  if t>0 
   then (re:sqrt(t/2), im:imagpart(z)/sqrt(2*t))
   else  (im:abs(z), re:0),
  return(float(re+im*%i))
 )$
 /*-----------------------------------*/ 

 Psi_n(r,t,z_last, Max_R):=
 /*   */
 block(
  [iMax:200,
  iMax2:0],
  /* -----  forward iteration of 2 points : z_last and w --------------*/
  array(forward,iMax-1), /* forward orbit of z_last for comparison */
  forward[0]:z_last,
  i:0,
  while cabs(forward[i])<Max_R  and  i< ( iMax-2) do
  (	
  /* forward iteration of z in fc plane & save it to forward array */
  forward[i+1]:forward[i]*forward[i] + c, /* z*z+c */
  /* forward iteration of w in f0 plane :  w(n+1):=wn^2 */
  r:r*2, /* square radius = R^2=2^(2*r) because R=2^r */
  t:mod(2*t,1),
  /* */
  iMax2:iMax2+1,
  i:i+1
  ),
  /* compute last w point ; it is equal to z-point */
  R:2^r,
  /* w:R*exp(2*%pi*%i*t),	z:w, */
  array(backward,iMax-1),
  backward[iMax2]:rectform(ev(R*exp(2*%pi*%i*t))), /* use last w as a starting point for backward iteration to new z */
  /* -----  backward iteration point  z=w in fc plane --------------*/
  for i:iMax2 step -1 thru 1 do
  (
  temp:csqrt(backward[i]-c), /* sqrt(z-c) */
  scalar_product:realpart(temp)*realpart(forward[i-1])+imagpart(temp)*imagpart(forward[i-1]),
  if (0>scalar_product) then temp:-temp, /* choose preimage */
  backward[i-1]:temp
  ),
  return(backward[0])
 )$

 GiveRay(t,c):=
 block(
  [r],
  /* range for drawing  R=2^r ; as r tends to 0 R tends to 1 */
  rMin:1E-10, /* 1E-4;  rMin > 0  ; if rMin=0 then program has infinity loop !!!!! */
  rMax:2, 
  caution:0.9330329915368074, /* r:r*caution ; it gives smaller r */
  /* upper limit for iteration */
  R_max:300,
  /* */
  zz:[], /* array for z points of ray in fc plane */
  /*  some w-points of external ray in f0 plane  */
  r:rMax,
  while 2^r<R_max do r:2*r, /* find point w on ray near infinity (R>=R_max) in f0 plane */
  R:2^r,
  w:rectform(ev(R*exp(2*%pi*%i*t))),
  z:w, /* near infinity z=w */
  zz:cons(z,zz),
  unless r<rMin do
  (	/* new smaller R */
  r:r*caution,  
  R:2^r,
  /* 
  w:rectform(ev(R*exp(2*%pi*%i*t))),*/
  /* */
  last_z:z,
  z:Psi_n(r,t,last_z,R_max), /* z=Psi_n(w) */
  zz:cons(z,zz)
  ),
  return(zz)
 )$

 /* Gives points of backward orbit of z=repellor       */
 GiveBackwardOrbit(c,repellor,zxMin,zxMax,zyMin,zyMax,iXmax,iYmax):=
  block(
   hit_limit:4, /* proportional to number of details and time of drawing */
   PixelWidth:(zxMax-zxMin)/iXmax,
   PixelHeight:(zyMax-zyMin)/iYmax,
   /* 2D array of hits pixels . Hit > 0 means that point was in orbit */
   array(Hits,fixnum,iXmax,iYmax), /* no hits for beginning */
  /* choose repeller z=repellor as a starting point */
  stack:[repellor], /*save repellor in stack */
  /* save first point to list of pixels  */ 
  x_y:[repellor], 
 /* reversed iteration of repellor */
  loop,
  /* pop = take one point from the stack */
  z:last(stack),
  stack:delete(z,stack),
  /*inverse iteration - first preimage (root) */
  z:finverseplus(z,c),
  /* translate from world to screen coordinate */
  iX:fix((realpart(z)-zxMin)/PixelWidth),
  iY:fix((imagpart(z)-zyMin)/PixelHeight),
  hit:Hits[iX,iY],
  if hit<hit_limit   
   then 
    (
    Hits[iX,iY]:hit+1,
    stack:endcons(z,stack), /* push = add z at the end of list stack */
    if hit=0 then x_y:endcons( z,x_y)
    ),
  /*inverse iteration - second preimage (root) */
  z:-z,
 /* translate from world to screen coordinate, coversion to integer */
  iX:fix((realpart(z)-zxMin)/PixelWidth),
  iY:fix((imagpart(z)-zyMin)/PixelHeight),
  hit:Hits[iX,iY],
  if hit<hit_limit   
   then 
    (
     Hits[iX,iY]:hit+1,
     stack:endcons(z,stack), /* push = add z at the end of list stack to continue iteration */
     if hit=0 then x_y:endcons( z,x_y)
    ),
   if is(not emptyp(stack)) then go(loop), 
 return(x_y) /* list of pixels in the form [z1,z2] */
 )$

compile(all);

 /* ----------------------- main ----------------------------------------------------*/

 start:elapsed_run_time ();
  c:-1; /* center of period 2 component */

 /* external angle in turns */
 /* resolution is proportional to number of details and time of drawing */
 iX_max:1000;
 iY_max:1000;
 /* define z-plane ( dynamical ) */
 ZxMin:-2.0;
 ZxMax:2.0;
 ZyMin:-2.0;
 ZyMax:2.0;

 /* compute ray points & save to zz list */
 zz1:GiveRay(1/3,c)$ 
 zz2:GiveRay(2/3,c)$
 
 /* limit cycle */
 z0:0;
 zp:[];
 zp:cons(z0,zp);
 z1:f(z0,c);
 zp:cons(z1,zp);
 z2:f(z1,c);
 zp:cons(z2,zp);

 /* compute fixed points */
 beta:rectform((1+csqrt(1-4*c))/2); /* compute repelling fixed point beta */
 alfa:rectform((1-csqrt(1-4*c))/2); /* other fixed point */

 /* compute backward orbit of repelling fixed point */
 xy: GiveBackwardOrbit(c,beta,ZxMin,ZxMax,ZyMin,ZyMax,iX_max,iY_max)$ /**/

 /* time of computations */
 time:fix(elapsed_run_time ()-start);

 /* draw it using draw package by */

 path:"~/maxima/"$ /* pwd, ended with / begining with ~ ,  if empty then file is in a home dir */

 load(draw); 
 draw2d(
  terminal  = png,
  file_name = sconcat(path,"J"),
  user_preamble="set size square;set key bottom right",
  title= concat("Dynamical plane for fc(z)=z*z+", string(c),"; Julia set and external rays landing on  fixed point z=alfa"),
  dimensions  = [iX_max, iY_max],
  yrange = [ZyMin,ZyMax],
  xrange = [ZxMin,ZyMax],
  xlabel     = "Z.re ",
  ylabel     = "Z.im",
  point_type = filled_circle,
  points_joined =true,
  point_size    =  0.2,
  color         = red,
  key = concat("external ray for angle ",string(1/3)),
  points(map(realpart,zz1),map(imagpart,zz1)),
  key = concat("external ray for angle ",string(2/3)),
  points(map(realpart,zz2),map(imagpart,zz2)),
  
  points_joined =false,
  color         = black,
  key = "backward orbit of z=beta",
  points(map(realpart,xy),map(imagpart,xy)),
  color         = blue,
  point_size    =  0.9,
  key = "repelling fixed point z= beta",
  points([[realpart(beta),imagpart(beta)]]),
  color         = yellow,
  key = "repelling fixed point z= alfa",
  points([[realpart(alfa),imagpart(alfa)]]),
  color         = green,
  key = "periodic z-points",
  points(map(realpart,zp),map(imagpart,zp))
 );