/*  
batch file for Maxima CAS
 */ 

start:elapsed_run_time ();

kill(all);
remvalue(all);

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

/* */
fn(p, z, c) :=
  if p=0 then z
  elseif p=1 then f(z,c)
  else f(fn(p-1, z, c),c);

/*Standard polynomial F_p \, which roots are periodic z-points of period p and its divisors */
F(p, z, c) := fn(p, z, c) - z ;

/* Function for computing reduced polynomial G_p\, which roots are periodic z-points of period p without its divisors*/
G[p,z,c]:=
block(
[f:divisors(p),
t:1], /* t is temporary variable = product of Gn for (divisors of p) other than p */
f:delete(p,f), /* delete p from list of divisors */
if p=1
then return(F(p,z,c)),
for i in f do 
 t:t*G[i,z,c],
g: F(p,z,c)/t,
return(ratsimp(g))
)$

GiveRoots(g):=
 block(
 [cc:bfallroots(expand(%i*g)=0)],
 cc:map(rhs,cc),/* remove string "c=" */
 cc:map('float,cc),
 return(cc)
  )$ 

/* 
circle D={w:abs(w)=1 } where w=l(t,r) 
t is angle in turns ; 1 turn = 360 degree = 2*Pi radians 
r is a radius 

*/
GiveC(angle,radius):=
(
 [w],  /* point of  unit circle   w:l(internalAngle,internalRadius); */
 w:radius*%e^(%i*angle*2*%pi),  /* point of circle */
 float(rectform(w/2-w*w/4))    /* point in a period 1 component of Mandelbrot set */
)$

/* endcons the complex point to list in the format for draw package */ 
endconsD(point,list):=endcons([realpart(point),imagpart(point)],list)$
consD(point,list):=cons([realpart(point),imagpart(point)],list)$

GiveForwardOrbit(z0,c,iMax):=
   /* 
   computes (without escape test)
   forward orbit of point z0
   and saves it to the list for draw package */
block(
 [z,orbit,temp],
 z:z0, /* first point = critical point z:0+0*%i */
 orbit:[[realpart(z),imagpart(z)]], 
 for i:1 thru iMax step 1 do
        ( z:expand(f(z,c)),
          orbit:endcons([realpart(z),imagpart(z)],orbit)),
         
 return(orbit) 
)$

/* gives 3 sublists from forward orbit of internal point */
GiveInternalRays(z0,c,iMax):= block
([a,b,d,z],
 a:[],
 b:[],
 d:[],
 z:z0,
 for i:1 thru iMax step 1 do
   ( 
   a:consD(z,a),
   z:f(z,c),
   b:consD(z,b),
   z:f(z,c),
   d:consD(z,d),
   z:f(z,c)
   ),
return([a,b,d])
)$

 /* 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] */
 )$

 
 
 /*-----------------------------------*/ 
 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:float(rectform(sqrt(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)
 )$

/* 

find symmetric point z3
 z3 is the same line as z1 and z2 such z2 is between z1 and z3
*/

GiveNextPoint(z1,z2):=(
[x,y,dx,dy],
 dx:realpart(z1)-realpart(z2),
 dy:imagpart(z1)-imagpart(z2),
 x:realpart(z2)-dx,
 y:imagpart(z2)-dy,
x+y*%i
)$

compile(all)$

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



period:3$

  

 /* 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;

 /* limit cycle */
 k:G[period,z,c]$ /* here c and z are symbols */

/* c:-0.12256+0.74486*%i;  value by Milnor*/
 c:0.74486176661974*%i-0.12256116687665; /* center of period 3 component */

 /* find periodic z points */
 s:GiveRoots(ev(k))$ /* ev moves value to c symbol here */ 
 z0:s[1];
 z1:rectform(float(f(z0,c)));
 z2:rectform(float(f(z1,c)));
 /* create 2 sublists : s1 and s2  from one list s */
 s1:[z0,z1,z2]$
 s2:delete(s[1],s);
 for z in s2 do if abs(z-z1)<0.1 then s2:delete(z,s2) ;
 for z in s2 do if abs(z-z2)<0.1 then s2:delete(z,s2) ;

 /* compute fixed points */
 beta:float(rectform((1+sqrt(1-4*c))/2)); /* compute repelling fixed point beta */
 alfa:float(rectform((1-sqrt(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)$ /**/


  /* compute ray points & save to zz list */
 eRay1o7:GiveRay(1/7,c)$
 eRay2o7:GiveRay(2/7,c)$
 eRay4o7:GiveRay(4/7,c)$  

 

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

 /* draw it using draw package by */
 load(draw); 

path:"~/maxima/batch/julia/rabbit/"$ /*  if empty then file is in a home dir */

/* if graphic  file is empty (= 0 bytes) then run draw2d command again */
 draw2d(
  terminal  = 'svg,
  file_name = sconcat(path,"Julia_1_3g"),
  user_preamble="set size square;set key bottom right",
  title= concat("Dynamical plane for fc(z)=z*z+",string(c)),
  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,
    
  
  points_joined =false,
  color         = black,
  key = "backward orbit of z=beta",
  points(map(realpart,xy),map(imagpart,xy)),
  

  points_joined =true,
  point_size    =  0.2,
  color         = red,
  key = "external ray 1/7",
  points(map(realpart,eRay1o7),map(imagpart,eRay1o7)),
  key = "external ray 2/7",
  points(map(realpart,eRay2o7),map(imagpart,eRay2o7)),
  key = "external ray 4/7",
  points(map(realpart,eRay4o7),map(imagpart,eRay4o7)),

  

  points_joined =false,
  
  color         = blue,
  point_size    =  1.4,
  key = "repelling fixed point z= beta",
  points([[realpart(beta),imagpart(beta)]]),
  color         = yellow,
  key = "fixed point alfa and repelling period 3 cycle",
  points([[realpart(alfa),imagpart(alfa)]]),
  color         = green,
  key = sconcat("attracting period ",string(period)," cycle"),
  points(map(realpart,s1),map(imagpart,s1))
 
 );