# Theory

Let $\Psi _{c}\,$  be the map ( Conformal map, isomorphism) from the complement (exterior) of the closed unit disk ${\overline {\mathbb {D} }}$  to the complement of the filled Julia set $\ K_{c}$ .

$\Psi _{c}:{\hat {\mathbb {C} }}\setminus {\overline {\mathbb {D} }}\to {\hat {\mathbb {C} }}\setminus K_{c}$

where ${\hat {\mathbb {C} }}$  denotes the Riemann sphere (extended complex plane.

Let $\Psi _{c}=\Phi _{c}^{-1}\,$  denote

The external ray of angle $\theta \,$  noted as ${\mathcal {R}}_{\theta }^{K}$ is:

• the image under $\Psi _{c}\,$  of straight lines ${\mathcal {R}}_{\theta }=\{\left(r\cdot e^{2\pi i\theta }\right):\ r>1\}$
${\mathcal {R}}_{\theta }^{K}=\Psi _{c}({\mathcal {R}}_{\theta })$
• set of points of exterior of filled-in Julia set with the same external angle $\theta$
${\mathcal {R}}_{\theta }^{K}=\{z\in {\hat {\mathbb {C} }}\setminus K_{c}:\arg(\Phi _{c}(z))=\theta \}$

# Important questions

## What are the unique rays?

• "The ray, t, which lands on the critical point (or on the root point of the bulb containing the critical point) of a Julia set is unique to that Julia set, so is a way of identifying the Julia set. So I call the julia set Jc by the name t. In the rest, I can only deal with julia sets where t=p/q is rational."
• Critically preperiodic polynomials are typically parameterized by the angle $\theta$  of the external ray landing at the critical value

## How to aproximate external rays?

In BDM images the external rays of angles (measured in turns):

$angle=(k/2^{n})~~{\mbox{mod }}~1\,$

can be seen as borders of subsets.

## when two external dynamic rays land at the same point

Questions:

• Q1: Given a point in the Julia set, how many rays are there landing at this point?
• Q2: Given two rays $R_{\theta _{1}}$  and $R_{\theta _{2}}$ , under what conditions do they land at the same point?
• When external rays land at the same point 
• The landing equivalence = Which rays land at the same point? 
• numerical criterion = they have the same ( numerically ) landing point
• Q3: Which rays support the same Fatou component?

Theorem 1.2 Let f be a polynomial of degree d ≥ 2 with locally connected Julia set. If two angles $\theta _{1}$  and $\theta _{2}$  have the same itinerary with respect to a critical portrait, then the landing points of external rays $R_{\theta _{1}}$  and $R_{\theta _{2}}$  either coincide or belong to the boundary of a Fatou domain, which is eventually iterated onto a Siegel disk.

## How to find angle of the external ray that land on the critical values z= c ?

That depends on the situation.

• There are combinatorial methods when the parameter c is a Misiurewicz point, for example.
• Most often you construct the angle first and then determine the parameter.
• But if you want to read an angle from the Julia set or from the Hubbard tree, you can proceed according to Demo 4 pages 1 & 2 from program Mandel

## Constructing the spine of filled Julia set

Algorithm for constructiong the spine is described by A. Douady

• join $-\beta \,$  and $\beta \,$ ,
• (to do )

## What means to draw external ray ?

It means:

• calculate (approximate) some points of ray
• join points by line segments

This will give an approximation of ray ${\mathcal {R}}$ .

# Drawing dynamic external ray

Algorithms

• Newton method
• backwards iteration
• inverse Boettcher map
• trace an external ray ( curve) for angle t in turns

## backwards iteration

This method has been used by several people and proved by Thierry Bousch.

Code in c++ by Wolf Jung can be found in procedure QmnPlane::backray() in file qmnplane.cpp ( see source code of the program Mandel ).

• Ray for periodic angle ( simplest case )

It will be explained by an example :

First choose external angle $t\,$  (in turns). External angle for periodic ray is a rational number.

$t={\frac {1}{3}}\,$

Compute period of external angle under doubling map.

Because "1/3 doubled gives 2/3 and 2/3 doubled gives 4/3, which is congruent to 1/3" 

$1\equiv 2{\pmod {3}}.\,$

or

${\frac {2}{3}}\equiv {\frac {1}{3}}{\pmod {1}}.\,$  

so external angle $t={\frac {1}{3}}\,$  has period 2 under doubling map.

on ray 1/3 is a point $w=R*e^{2\pi i/3}\,$

on ray 2/3 is a point $w=R*e^{4\pi i/3}\,$ .

Near infinity $z=w\,$  so one can swith to dynamical plane ( Boettcher conjugation )

Backward iteration (with proper chose from two possibilities) of point on ray 1/3 goes to ray 2/3, back to 1/3 and so on.

In C it is :

/* choose one of 2 roots: zNm1 or -zNm1  where zN = sqrt(zN - c )  */
if (creal(zNm1)*creal(zN) + cimag(zNm1)*cimag(zN) <= 0) zNm1=-zNm1;


or in Maxima CAS :

if (z1m1.z01>0) then z11:z1m1 else z11:-z1m1;


One has to divide set of points into 2 subsets ( 2 rays). Draw one of these 2 sets and join the points. It will be an approximation of ray.

• Ray for preperiodic angle ( to do )
/*

compute last point ~ landing point
of the dynamic ray for periodic angles ( in turns )

gcc r.c -lm -Wall -march=native

landing point of ray for angle = 1 / 15 = 0.0666666666666667 is = (0.0346251977103306 ;  0.4580500411138030 ) ; iDistnace = 18
landing point of ray for angle = 2 / 15 = 0.1333333333333333 is = (0.0413880816505388 ;  0.5317194187688231 ) ; iDistnace = 17
landing point of ray for angle = 4 / 15 = 0.2666666666666667 is = (-0.0310118081927549 ;  0.5440125864026020 ) ; iDistnace = 17
landing point of ray for angle = 8 / 15 = 0.5333333333333333 is = (-0.0449867688014234 ;  0.4662592852362425 ) ; iDistnace = 18

*/

#include <stdio.h>
#include <stdlib.h> // malloc
#include <math.h> // M_PI; needs -lm also
#include <complex.h>

/* --------------------------------- global variables and consts ------------------------------------------------------------ */
#define iPeriodChild 4 // iPeriodChild of secondary component joined by root point

// - --------------------- functions ------------------------

/*
principal square  root of complex number
wikipedia  Square_root

z1= I;
z2 = root(z1);
printf("zx  = %f \n", creal(z2));
printf("zy  = %f \n", cimag(z2));
*/
double complex root(double x, double y)
{

double u;
double v;
double r = sqrt(x*x + y*y);

v = sqrt(0.5*(r - x));
if (y < 0) v = -v;
u = sqrt(0.5*(r + x));
return u + v*I;
}

// This function only works for periodic  angles.
// You must know the iPeriodChild n before calling this function.
// draws all "iPeriodChild" external rays
// commons  File:Backward_Iteration.svg
// based on the code by Wolf Jung from program Mandel
// http://www.mndynamics.com/

int ComputeRays( //unsigned char A[],
int n, //iPeriodChild of ray's angle under doubling map
int iterMax,
double Cx,
double Cy,
double dAlfaX,
double dAlfaY,
double PixelWidth,
complex double zz[iPeriodChild] // output array

)
{
double xNew; // new point of the ray
double yNew;

const double R = 10000; // very big radius = near infinity
int j; // number of ray
int iter; // index of backward iteration

double t,t0; // external angle in turns
double num, den; // t = num / den

double complex zPrev;
double u,v; // zPrev = u+v*I

int iDistance ; // dDistance/PixelWidth = distance to fixed in pixels

/* dynamic 1D arrays for coordinates ( x, y) of points with the same R on preperiodic and periodic rays  */
double *RayXs, *RayYs;
int iLength = n+2; // length of arrays ?? why +2

//  creates arrays :  RayXs and RayYs  and checks if it was done
RayXs = malloc( iLength * sizeof(double) );
RayYs = malloc( iLength * sizeof(double) );
if (RayXs == NULL || RayYs==NULL)
{
fprintf(stderr,"Could not allocate memory");
getchar();
return 1; // error
}

// external angle of the first ray
num = 1.0;
den = pow(2.0,n) -1.0;
t0 = num/den; // http://fraktal.republika.pl/mset_external_ray_m.html
t=t0;
// printf(" angle t = %.0f / %.0f = %f in turns \n", num, den, t0);

//  starting points on preperiodic and periodic rays
//  with angles t, 2t, 4t...  and the same radius R
for (j = 0; j < n; j++)
{ // z= R*exp(2*Pi*t)
RayXs[j] = R*cos((2*M_PI)*t);
RayYs[j] = R*sin((2*M_PI)*t);
//
// printf(" %d angle t = = %.0f / %.0f =  %.16f in turns \n", j, num , den,  t);
//
num *= 2.0;
t *= 2.0; // t = 2*t
if (t > 1.0) t--; // t = t modulo 1

}
//zNext = RayXs + RayYs *I;

// printf("RayXs  = %f \n", RayXs);
// printf("RayYs  = %f \n", RayYs);

// z[k] is n-periodic. So it can be defined here explicitly as well.
RayXs[n] = RayXs;
RayYs[n] = RayYs;

//   backward iteration of each point z
for (iter = -10; iter <= iterMax; iter++)
{

for (j = 0; j < n; j++) // n +preperiod
{ // u+v*i = sqrt(z-c)   backward iteration in fc plane
zPrev = root(RayXs[j+1] - Cx , RayYs[j+1] - Cy ); // , u, v
u=creal(zPrev);
v=cimag(zPrev);

// choose one of 2 roots: u+v*i or -u-v*i
if (u*RayXs[j] + v*RayYs[j] > 0)
{ xNew = u; yNew = v; } // u+v*i
else { xNew = -u; yNew = -v; } // -u-v*i

// draw part of the ray = line from zPrev to zNew
// dDrawLine(A, RayXs[j], RayYs[j], xNew, yNew, j, 255);

//
RayXs[j] = xNew; RayYs[j] = yNew;

} // for j ...

//RayYs[n+k] cannot be constructed as a preimage of RayYs[n+k+1]
RayXs[n] = RayXs;
RayYs[n] = RayYs;

// convert to pixel coordinates
//  if z  is in window then draw a line from (I,K) to (u,v) = part of ray

// printf("for iter = %d cabs(z) = %f \n", iter, cabs(RayXs + RayYs*I));

}

// Approximate end of ray by straight line to it's landing point here = alfa fixed point
// for (j = 0; j < n + 1; j++)
//  dDrawLine(A, RayXs[j],RayYs[j], dAlfaX, dAlfaY,j, 255 );

// this check can be done only from inside this function
t=t0;
num = 1.0;
for (j = 0; j < n ; j++)
{

zz[j] = RayXs[j] +  RayYs[j] * I; // save to the output array
// Approximate end of ray by straight line to it's landing point here = alfa fixed point
//dDrawLine(RayXs[j],RayYs[j], creal(alfa), cimag(alfa), 0, data);
iDistance = (int) round(sqrt((RayXs[j]-dAlfaX)*(RayXs[j]-dAlfaX) +  (RayYs[j]-dAlfaY)*(RayYs[j]-dAlfaY))/PixelWidth);
printf("last point of the ray for angle = %.0f / %.0f = %.16f is = (%.16f ;  %.16f ) ; Distance to fixed = %d  pixels \n",num, den, t, RayXs[j], RayYs[j],  iDistance);
num *= 2.0;
t *= 2.0; // t = 2*t
if (t > 1) t--; // t = t modulo 1
} // end of the check

// free memmory
free(RayXs);
free(RayYs);

return  0; //
}

int main()
{

complex double l[iPeriodChild];
int i;
// external angle in turns = num/den;
double num = 1.0;
double den = pow(2.0, iPeriodChild) -1.0;

ComputeRays( iPeriodChild,
10000,
0.25, 0.5,
0.00, 0.5,
0.003,
l ) ;

printf("\n see what is in the output array : \n");

for (i = 0; i < iPeriodChild ; i++) {
printf("last point of the ray for angle = %.0f / %.0f = %.16f is = (%.16f ;  %.16f ) \n",num, den, num/den, creal(l[i]), cimag(l[i]));
num *= 2.0;}

return 0;
}


Run it :

./a.out


And output :

last point of the ray for angle = 1 / 15 = 0.0666666666666667 is = (0.0346251977103306 ;  0.4580500411138030 ) ; Distance to fixed = 18  pixels
last point of the ray for angle = 2 / 15 = 0.1333333333333333 is = (0.0413880816505388 ;  0.5317194187688231 ) ; Distance to fixed = 17  pixels
last point of the ray for angle = 4 / 15 = 0.2666666666666667 is = (-0.0310118081927549 ;  0.5440125864026020 ) ; Distance to fixed = 18  pixels
last point of the ray for angle = 8 / 15 = 0.5333333333333333 is = (-0.0449867688014234 ;  0.4662592852362425 ) ; Distance to fixed = 19  pixels

see what is in the output array :
last point of the ray for angle = 1 / 15 = 0.0666666666666667 is = (0.0346251977103306 ;  0.4580500411138030 )
last point of the ray for angle = 2 / 15 = 0.1333333333333333 is = (0.0413880816505388 ;  0.5317194187688231 )
last point of the ray for angle = 4 / 15 = 0.2666666666666667 is = (-0.0310118081927549 ;  0.5440125864026020 )
last point of the ray for angle = 8 / 15 = 0.5333333333333333 is = (-0.0449867688014234 ;  0.4662592852362425 )



Point on the ray is moving backwards:

• very fast when it it far from Julia set
• very slow near Julia set ( after 50 iterations distance between points = 0 pixels )

See example computing ( here pixel size = 0.003 ) :

# iteration distance_between_points_in_pixels
0 	 3300001
1 	 30007
2 	 2296
3 	 487
4 	 179
5 	 92
6 	 54
7 	 34
8 	 23
9 	 18
10 	 14
11 	 11
12 	 9
13 	 7
14 	 6
15 	 5
16 	 5
17 	 4
18 	 4
19 	 3
20 	 3
21 	 3
22 	 3
23 	 2
24 	 2
25 	 2
26 	 2
27 	 2
28 	 2
29 	 2
30 	 1
31 	 1
32 	 1
33 	 1
34 	 1
35 	 1
36 	 1
37 	 1
38 	 1
39 	 1
40 	 1
41 	 1
42 	 1
43 	 1
44 	 1
45 	 1
46 	 1
47 	 1
48 	 1
49 	 1
50 	 1
51 	 1
52 	 1
53 	 1
54 	 1
55 	 1
56 	 0
57 	 0
58 	 0
59 	 0
60 	 0


One can choose points which differ by pixel size :

#iteration    distance(z1,z2)   distance (z,alfa)
0 	 3300001 	 33368
1 	   30007 	  3364
2 	    2296 	  1074
3 	     487 	   591
4 	     179 	   413
5 	      92 	   321
6 	      54 	   267
7 	      34 	   234
8 	      23 	   211
9 	      18 	   193
10 	      14 	   179
11 	      11 	   169
12 	       9 	   160
13 	       7 	   153
14 	       6 	   146
15 	       5 	   141
16 	       5 	   136
17 	       4 	   132
18 	       4 	   128
19 	       3 	   125
20 	       3 	   122
21 	       3 	   119
22 	       3 	   117
23 	       2 	   115
24 	       2 	   112
25 	       2 	   110
26 	       2 	   109
27 	       2 	   107
28 	       2 	   105
29 	       2 	   104
30 	       1 	   102
31 	       1 	   101
32 	       1 	   100
33 	       1 	    99
34 	       1 	    97
35 	       1 	    96
36 	       1 	    95
38 	       2 	    93
40 	       2 	    92
42 	       2 	    90
44 	       2 	    88
46 	       1 	    87
48 	       1 	    86
50 	       1 	    84
52 	       1 	    83
54 	       1 	    82
56 	       1 	    81
59 	       1 	    80
62 	       1 	    78
65 	       1 	    77
68 	       1 	    76
71 	       1 	    75
74 	       1 	    74
78 	       1 	    73
82 	       1 	    71
86 	       1 	    70
90 	       1 	    69
95 	       1 	    68
100 	       1 	    67
105 	       1 	    66
111 	       1 	    65
117 	       1 	    64
124 	       1 	    63
131 	       1 	    62
139 	       1 	    60
147 	       1 	    59
156 	       1 	    58
166 	       1 	    57
177 	       1 	    56
189 	       1 	    55
202 	       1 	    54
216 	       1 	    53
231 	       1 	    52
247 	       1 	    51
265 	       1 	    50
285 	       1 	    49
307 	       1 	    48
331 	       1 	    47
358 	       1 	    46
388 	       1 	    45
421 	       1 	    44
458 	       1 	    43
499 	       1 	    42
545 	       1 	    41
597 	       1 	    40
655 	       1 	    39
721 	       1 	    38
796 	       1 	    37
881 	       1 	    36
978 	       1 	    35
1090 	       1 	    34
1219 	       1 	    33
1368 	       1 	    32
1542 	       1 	    31
1746 	       1 	    30
1986 	       1 	    29
2270 	       1 	    28
2608 	       1 	    27
3013 	       1 	    26
3502 	       1 	    25
4098 	       1 	    24
4830 	       1 	    23
5737 	       1 	    22
6873 	       1 	    21
8312 	       1 	    20
10157 	       1 	    19
12555 	       1 	    18
15719 	       1 	    17
19967 	       1 	    16
25780 	       1 	    15
33911 	       1 	    14
45574 	       1 	    13
62798 	       1 	    12
89119 	       1 	    11
131011 	       1 	    10
201051 	       1 	     9
325498 	       1 	     8
564342 	       1 	     7
1071481 	       1 	     6
2308074 	       1 	     5
5996970 	       1 	     4
21202243 	       1 	     3
136998728 	       1 	     2


One can see that moving from pixel 12 to 11 near alfa needs 27 000 iterations. Computing points up to 1 pixel near alfa needs : 2m1.236s

## Drawing dynamic external ray using inverse Boettcher map by Curtis McMullen

This method is based on C program by Curtis McMullen and its Pascal version by Matjaz Erat

There are 2 planes here :

• w-plane ( or f0 plane )
• z-plane ( dynamic plane of fc plane )

Method consist of 3 big steps :

• compute some w-points of external ray of circle for angle $\theta \,$  and various radii (rasterisation)
$w_{0}=R*e^{2\pi it}\,$  where $1
• map w-points to z-point using inverse Boettcher map $\Psi _{c}=\Phi _{c}^{-1}\,$
$z_{0}=\Psi _{c}(w_{0})=f_{c}^{-n}(w_{0}^{2^{n}})\,$
• draw z-points ( and connect them using segments ( line segment is a part of a line that is bounded by two distinct end points )

First and last steps are easy, but second is not so needs more explanation.

### Rasterisation

For given external ray in $f_{0}\,$  plane each point of ray has :

• constant value $t\,$  ( external angle in turns )
• variable radius $R\,$

so $w\,$  points of ray are parametrised by radius $R\,$  and can be computed using exponential form of complex numbers :

$w=F_{t}(R)=R*e^{2*\Pi *i*t}\,$

One can go along ray using linear scale :

t:1/3; /* example value */
R_Max:4;
R_Min:1.1;
for R:R_Max step -0.5 thru R_Min do w:R*exp(2*%pi*%i*t);
/* Maxima allows non-integer values in for statement */


It gives some w points with equal distance between them.

Another method is to use nonlinera scale.

To do it we introduce floating point exponent $r\,$  such that :

$R=2^{r}\,$

and

$w=G_{t}(r)=2^{r}*e^{2*\Pi *i*t}\,$

To compute some w points of external ray in $f_{0}\,$  plane for angle $t\,$  use such Maxima code :

t:1/3; /* external angle in turns */
/* range for computing  R ; as r tends to 0 R tends to 1 */
rMax:2; /* so Rmax=2^2=4 /
rMin:0.1; /* rMin > 0   */
caution:0.93; /* positive number < 1 ; r:r*caution  gives smaller r */
r:rMax;
unless r<rMin do
(
r:r*caution, /* new smaller r */
R:2^r, /* new smaller R */
w:R*exp(2*%pi*%i*t) /* new point w in f0 plane */
);


In this method distance between points is not equal but inversely proportional to distance to boundary of filled Julia set.

It is good because here ray has greater curvature so curve will be more smooth.

### Mapping

• forward iteration in $f_{0}\,$  plane
$w_{n}=w_{0}^{2^{n}}=R^{2^{n}}*e^{2\pi i2^{n}t}\,$

until $w_{n}\,$  is near infinity

$abs(w_{n})>R_{limit}\,$
• switching plane ( from $f_{0}\,$  to $f_{c}\,$ )
$z_{n}=w_{n}\,$

( because here, near infinity : $w=z\,$ )

• backward iteration in $f_{c}\,$  plane the same ( $n\,$ ) number of iterations
• last point $z_{0}\,$  is on our external ray

1,2 and 4 minor steps are easy. Third is not.

$z_{n-1}=f_{c}^{-1}(z_{n})={\sqrt {z_{n}-c}}\,$

Backward iteration uses square root of complex number. It is 2-valued functions so backward iteration gives binary tree.

One can't choose good path in such tree without extre informations. To solve it we will use 2 things :

• equicontinuity of basin of attraction of infinity
• conjugacy between $f_{0}\,$  and $f_{c}\,$  planes

### Equicontinuity of basin of attraction of infinity

Basin of attraction of infinity ( complement of filled-in Julia set) contains all points which tends to infinity under forward iteration.

$A_{f_{c}}(\infty )\ {\overset {\underset {\mathrm {def} }{}}{=}}\ \{z\in \mathbb {C} :f_{c}^{(k)}(z)\to \infty \ as\ k\to \infty \}.$

Infinity is superattracting fixed point and orbits of all points have similar behaviour. In other words orbits of 2 points are assumed to stay close if they are close at the beginning.

It is equicontinuity ( compare with normality).

In $f_{c}\,$  plane one can use forward orbit of previous point of ray for computing backward orbit of next point.

### Detailed version of algorithm

• compute first point of ray (start near infinity ang go toward Julia set )
$w=2^{r}*e^{2\Pi it}\,$  where $r=r_{max}\,$

here one can easily switch planes :

$z=w\,$

It is our first z-point of ray.

• compute next z-point of ray
• compute next w-point of ray for $r=r*caution\,$
• compute forward iteration of 2 points : previous z-point and actual w-point. Save z-orbit and last w-point
• switch planes and use last w-point as a starting point :$z_{n}=w_{last}\,$
• backward iteration of new $z_{n}\,$  toward new $z_{0}\,$  using forward orbit of previous z point
• $z_{0}\,$  is our next z point of our ray
• and so on ( next points ) until $r

Maxima CAS src code

/* gives a list of z-points of external ray for angle t in turns and coefficient c */
GiveRay(t,c):=
block(
[r],
/* range for drawing  R=2^r ; as r tends to 0 R tends to 1 */
rMin:1E-20, /* 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)
)\$