# Fractals/Iterations in the complex plane/MandelbrotSetExterior/ParameterExternalRay

Parameter external ray
Tuning external ray

## landing

The Douady-Hubbard landing theorem for periodic external for a complex polynomial f with bounded postcritical set:
* every periodic external ray lands at a repelling or parabolic periodic point
* conversely every repelling or parabolic point is the landing point of at least one periodic external ray

# Q&A
## relation between parameter and dynamic rays
### wake

Mandel demo 3 ( external ray) page 8
For a p-periodic angle, the corresponding dynamic ray of a connected Julia set lands at a repelling or parabolic point with period dividing p. The corresponding parameter ray lands at the root of a hyperbolic component with period p.
To show this, consider a parameter c0 where the parameter ray accumulates. For parameters c on this ray, the corresponding dynamic ray branches. If the dynamic ray for c0 landed at a repelling periodic point, it would land without branching also for c ≈ c0 by the Implicit Function Theorem and a pullback argument. This contradiction shows that the dynamic ray for the parameter c0 must land at a parabolic point, and c0 is a root.
In fact, exactly two parameter rays with periodic angles land at every root. This was shown
* by Adrien Douady and John Hamal Hubbard, using Fatou coordinates for the parabolic implosion (Chapter 2).
* Combinatorial proofs were given by Dierk Schleicher and John Milnor.
For parameters c in the wake between the two parameter rays:
* the two dynamic rays with the same angles land together,
* and the critical value z = c is always in the wake between them, even if there are more dynamic rays landing.
The left image shows the parameter rays with the 2-periodic angles 1/3 and 2/3, which land at the root of the period-2 component. The cross indicates the center c = -1.
The right image shows the corresponding filled Julia set, the "basilica," and the two dynamic rays landing at the repelling fixed point αc . The cross is at the critical value z = c. Now you shall see the Julia sets for the centers corresponding to the p-periodic angles 1/(2p-1) and 2/(2p-1), and all parameter rays with these denomiantors: Hit Return or push the Go-button to increase the period (2...9).

### Subdivision rule

How to compute dynamic external rays that land on the critical point $z=0$ of the dendrite Julia set?
Steps:
* the preperiodic angle $\theta$ (decimal fraction with even denominator )
* find on the parameter plane the landing point of the parameter external ray $\theta$ which is Misiurewicz point: $\gamma _{M}(\theta )=c=M$
* on the dynamic plane:
** landing point of the ray with angle $\theta$ is critical value: $z=c$
** on the critical point z = 0 land rays with angles: $\theta /2$ and $(\theta +1)/2$
Examples of parameter rays:
* Ray for angle ${\frac {1}{2}}=0.0(1)$ lands on the point $c=\gamma _{M}(1/2)=-2=M_{1,1}$ from the parameter plane. It is tip of the main antenna ( end of 1/2 limb).
* Ray for angle ${\frac {1}{4}}=0.00(1)$ lands on the point $c=\gamma _{M}(1/4)=-0.228155493653962+1.115142508039937i=M_{2,1}$ from the parameter plane. It is the first tip of wake 1/3.
* Ray for angle ${\frac {1}{6}}=0.0(01)$ lands on the point $c=\gamma _{M}(1/6)=i=M_{1,2}$ from the parameter plane. It is the last tip of wake 1/3
* Ray for angle ${\frac {17}{240}}=0001(0010)$ lands on the point $c=\gamma _{M}(17/240)=0.366362983422764+0.591533773261445i=M_{4,4}$ from the parameter plane. It is the principle Misiurewicz point ( branch point or hub) of wake 1/4.
* Ray for angle ${\frac {9}{56}}=0.001(010)$ lands on the point $c=\gamma _{M}(9/56)=-0.101096363845622+0.956286510809142i=M_{3,3}$ from the parameter plane. It is the principle Misiurewicz point ( branch point or hub) of wake 1/3.
* Ray for angle ${\frac {129}{16256}}=0.0000001(0000010)$ lands on the point $c=\gamma _{M}(129/16256)=0.397391822296541+0.133511204871878i=M_{7,7}$ from the parameter plane. It is the principle Misiurewicz point ( branch point or hub) of wake 1/7.
Partition of the dynamic plane by dynamic rays is related with the kneading sequence
* 1/4
* 1/6
* 9/56
* 129/16256

### principal Misiurewicz point rule

Notes from demo 3 (external rays) page 9/12 from program mandel by Wolf Jung
The parameter rays with the angles 1/7 and 2/7 land at the root of a period-3 component, which is of satellite type with rotation number 1/3.
For all parameters c in the wake between the rays
* for c in the 1/3-limb of the Mandelbrot set,
** for c equal to the principle Misiurewicz point of the wake
the dynamic rays with the angles 1/7, 2/7, and 4/7:
* land together at the repelling fixed point $z=\alpha _{c}$
* the critical value z = c is between the first two rays.
We shall compute the external angles of certain preimages of fixed point alpha $z=\alpha _{c}$under $f_{c}(z)$.
Note that an angle $\theta$ has two preimages under doubling modulo 1, $(\theta /2,(\theta +1)/2)$.
Point $z=-\alpha _{c}$ is the only preimage of $z=\alpha _{c}$ different from the fixed point itself.
The angle 1/7 has the preimages (1/7)/2 = 1/14 and (1/7 + 1)/2 = 4/7. The latter angle belongs to $z=\alpha _{c}$, so 1/14 is an external angle of $z=-\alpha _{c}$ .
In the same way, the other angles 9/14 and 11/14 are obtained. The rays are drawn blue.
Move z to that preimage of $z=-\alpha _{c}$ between the rays for 2/7 and 4/7. The angle 1/14 has the preimages (1/14)/2 = 1/28 and (1/14 + 1)/2 = 15/28. Only the latter angle is in the chosen interval. The other two external angles of z are 9/28 and 11/28. The rays are drawn magenta.
Now z is the preimage between the rays for 1/7 and 2/7. By taking preimages in this interval, the external angles 9/56, 11/56, and 15/56 are obtained. The rays are drawn red.
Rays with preperiodic angles, i.e., even denominators, land at preperiodic points in the dynamic plane, or at Misiurewicz points in the parameter plane. For these parameters, the critical value is preperiodic under the iteration of fc(z).
orbit of repelling fixed point alpha $z=\alpha _{c}$ has period 3 :
| z | orbit portrait | description |
| -0.327586178699357 +0.577756453298949 i | (1/7, 2/7, 4/7) | repelling fixed point alpha $z=\alpha _{c}$ |
| +0.327586178699356 -0.577756453298949 i | (1/14, 9/14, 11/14) | preimage: $z=b(\alpha _{c})=f^{-1}(\alpha _{c})=-\alpha _{c}$ |
| -1.005359779818338 +0.762932332734299 i | (9/28, 11/28, 15/28 ) | a(z) |
| -0.101096363845622 +0.956286510809142 i | (9/56, 11/56, and 15/56) | a(z) |
| -0.000000003174826 +0.000000002063599 i | (9/56, 11/56, and 15/56) | a(z) ; it should be critical point z = 0 |
| -0.728941649258717 +0.655941740822478 i | | a(z) |
| -0.184581740009785 +0.813582020547489 i | | a(z) |
| -0.202294024682997 +0.352715534938011 i | | a(z) |
| -0.505370243253601 +0.597157217578291 i | | a(z) |
| -0.261224737974777 +0.687395259758146 i | | a(z) |
| -0.276433581450372 +0.486357789166200 i | | a(z) |

## symbolic dynamic

* Symbolic Dynamics of Quadratic Polynomials. Version of July 27, 2011. Henk Bruin, Alexandra Kaffl and Dierk Schleicher
** long version
** short version

## idea

* take a segment of straight ray ( near infinity )
* pull it back toward boundary of Mandelbrot set

## How to aproximate external rays?

* binary decomposition method for Julia set ( BDM/J)
* field lines of the scalar field ( potential)
* SAM = Stripe Average Method is basically a cheap way to calculate an angle.
* BDM, some rays can be seen
* boundaries of BDM
* field lines
* SAM
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.

## What means to draw external ray ?

It means:
* calculate (approximate) DS points of ray. The result is the set of complex numbers $\{c_{m}:1\leq m\leq DS\}$ ( points on the parameter plane ), use numerical algotrithm
* join points by line segments, use graphical algorithm )
This will give an approximation of ray ${\mathcal {R}}$.
trace an external ray for angle t in turns, which means ( description by Claude Heiland-Allen)
* starting at a large circle (e.g. r = 65536, x = r cos(2 pi t), y = r sin(2 pi t))
* following the escape lines (like the edges of binary decomposition colouring with large escape radius) in towards the Mandelbrot set.
this algorithm is O(period^2), which means doubling the period takes 4x as long, so it's only really feasible for periods up to a few 1000.

## How to compute one point of the ray ?

By solving polynomial equation
with numerical methods. The root of above equation is point $c_{m}$.
It is a point of the external parameter ray ${\mathcal {R}}(\theta )$
or
Using Newton method ( iteration ) one can compute approximation of point $c_{m}$
What one needs to start :
* arbitrarily chosen external angle $\theta$ of the ray $R(\theta )$ one wants to draw. Angle is usually given in turns
* value of function P ( which approximates Boettcher mapping ) and its derivative P'
* starting point $c_{m,0}$ ( initial approximation )
* stopping rule ( criteria to stop iteration ): Ray tracing has a natural stopping condition: when the ray enters the atom domain with period p, Newton's method is very likely to converge to the nucleus at its center.

## When ray lands ?

"The rays get closer and closer to the boundary, but don't reach it in finite time - for a more exact boundary point you need to switch to different methods when the ray is close enough. For points on the boundary of hyperbolic components, split the internal angled address (computable from the angle) into island and child path components, when tracing the ray to the parent island use atom domain test (to see if Newton's method is likely to converge to the right place) and switch to Newton's method to find the nucleus of the parent island and then trace internal rays through the chain of connected components to the desired boundary point.
For Misiurewicz points, there is probably a similar test to the atom domain test after which point Newton's method will converge to the desired location (though rays to Misiurewicz points tend to converge much more quickly than rays to roots of hyperbolic components anyway). The atom domain test checks that the iteration count of the last minimum of |z| is the same as the period of the ray." Claude Heiland-Allen
So ray does not "land" in the finite time. Landing point can be denoted as $c_{\infty }$

# Algorithms

* Jungreis-Ewing-Schober (JES) algorithm (based on the Laurent series expansion about infinity for $\Phi _{M}^{-1}\,$)
* OTIS method ( based on the Newton method)

# tracing rays

* outwards: "External Rays of the form 2pi*n/32, on top of the modulus of the potential gradient. For each point c, a path is created that follows the direction of the gradient of the potential. Each step size is proportional to the distance estimation to M. When the point is far enough of M, it's phase aproximates the phase of phi(c)." Inigo Quilez
* inwards : "The drawing method : ... the path is followed in reverse order: from the infinity towards M, following the minus gradient."
Tracing ray
* inwards = from infinity towards Mandelbrot set
* outwards = from point near Mandelbrot set towards infinity
Collecting bits of external angle's binary expansion when crossing dwell bands ( boundaries of level sets):
* inwards: add bit at the end of binary expansion
* outwards: add the bit at the beginning of binary expansion
crossing dwell bands ( boundaries of level sets)
* when ray value is changing it's integer part
See also
* binary decomposition
** parameter plane
** dynamic plane
* reading bits from binary decomposition = collecting bits of external angle's binary expansion

## in

* When tracing inwards, one peels off the most-significant bit (aka angle doubling) each time the ray crosses a dwell band (integer part of normalized iteration count increases by 1).
trace an external ray for angle t in turns, which means ( description by Claude Heiland-Allen)
* starting at a large circle (e.g. r = 65536, x = r cos(2 pi t), y = r sin(2 pi t))
* following the escape lines (like the edges of binary decomposition colouring with large escape radius) in towards the Mandelbrot set.
this algorithm is O(period^2), which means doubling the period takes 4x as long, so it's only really feasible for periods up to a few 1000.

## out
# Newton method

* using Newton method
* description by Tomoki Kawahira
* tracing inward ( from infinity toward Mandelbrot set) = ray-in
* arbitrary precision ( mpfr) with dynamic precision adjustment by Claude Heiland-Allen

## variables

* r = radial parameter = radius ( see complex potential )
* m = radial index = index of point along ray, integer
* j = sharpness = number of points on the dwell band, integer
* k = integer depth = number of dwell bands, integer
* d = m/S = real depth, floating point number ( name d is not used by Kawahira)
* l = index of Newton iteration ( name l is not used by Kawahira)
* n = index of iteration for computing Newton map
Names are from T Kawahira description

### Radial parameter r
## constant values

* S = $j_{max}$ =
* D = $k_{max}$ =
* R = ER = Escape Radius
* DS = $m_{max}$ = number of points
* $L=L_{m}=l_{m,max}$

## ranges

* sharpness : $1\leq j\leq S$
* radial index : $1\leq m\leq DS$
* depth :
** $1\leq k\leq D$
** $1\leq d\leq D$
* ray's radius ( subset): $R^{1/2^{D}}\leq r<R$
* iterations :
** quadratic : $1\leq n\leq k$
** Newton : $1\leq l_{m}\leq L_{m}$

## sequences

m-sequences ( along the ray toward the Mandelbrot set):
* $r_{1},r_{2},\cdots ,r_{SD}$
* $t_{1},t_{2},\cdots ,t_{SD}$
* $c_{1},c_{2},\cdots ,c_{SD}$
Newton sequences = l-sequences, here m is constant:
* from initial value $c_{m,0}$ toward an approximation of $c_{m,L}$
* $c_{m,0},c_{m,1},\cdots ,c_{m,L}$

## Maps
### r map

compare it with inverse iteration on c=0 dynamic plane

#### Depth

Using fixed integer D (maximal depth ) :
$k_{max}=D>1$
and fixed maximal value of radial parameter ( escape radius = ER ) :
$r_{max}=R>1$
one can compute D points of ray using formula :
$r=R^{1/2^{k}}$
which is :
$1<R^{1/2^{D}}<=r<=R$
When $k\to k_{max}=D\$ then $r\to 1$ and radius reaches enough close to the boundary of Mandelbrot set
Output :

#### Depth and sharpness

How to make ray more smooth ? Add more points between level sets.
Using:
* fixed integer S =maximal sharpness
* fixed integer D = maximal depth
one can compute S*D points of ray using fomula :
Note that k is equal to integer part of d :
and last point is the same as in depth method
but there are more points here because :
Output :

#### m

One can use only one loop : m-loop and ccompute j,k and d from m
output :

### Polynomial map q

Complex quadratic polynomial :
iteration :

### Map t

compare it with forward iteration on c=0 plane :

### Polynomial map P

P is a polynomial of degree $2^{k}$ in variable c.
Derivative with respect to c :

### Newton map N

Newton map:
note that here :

#### How to compute new value $c_{m,l+1}$ ?

* start with $c_{m,l}$
* compute value $t_{m}$ ( without iteration )
* compute $C_{n}$ and $D_{n}$ using quadratic iteration from n=1 to n=k
* compute next point $c_{m,l+1}$ using one Newton iteration
Arbitrary names :
Note that the derivative of a constant is zero.
The recursive formulae and initial values:
After k quadratic iterations compute new value $c_{m,l+1}$ using one Newton iteration
It is implemented in :
* mandelbrot-numerics library
** c/lib/m_d_exray_in_step = double precision
** c/lib/m_r_exray_in.c = mpfr ( arbitrary precision)

#### Newton iteration

Formula :
Newton iterations gives Newton sequence ( = l-sequence, here m is constant):
* from initial value $c_{m,0}$ toward an approximation of $c_{m,L}$
* $c_{m,0},c_{m,1},\cdots ,c_{m,L}$
Sequence :

#### Initial points

* $c_{1,0}=Re^{2\pi i\theta }$
* The value $c_{m-1,L}$ is presumably a “neighbor” of $c_{m,L}$ on ray so use it as the initial value for $c_{m}$ which is $c_{m,0}$

# Code

* c code by Claude Heiland-Allen
** book console program and gui version, uses stream
** mandelbrot-numerics library
** mandelbrot perturbator GUI program
* newtonRay - cpp code from program Mandel by Wolf Jung
* sage
** description by Ben Barros
** mandel_julia.py = python code at git repo by Ben Barros
** mandel_julia_helper.pyx
** final-report
** Sage Reference Manual: Discrete dynamics Release 8.3
* The Mandelbrot capacitor with field lines by Nils Berglund and c code : heat.c

## c code
### mpfr

Code for computing external ray inwards using Newton method based on the mandelbrot-book code by Claude Heiland-Allen
External angle is read from the string ( binary fraction)

## Python

Ray in procedure

# Test and precision

| angle $\theta$ in turns | landing point $c_{\infty }$ |
| 0 | 1/4 |
| 1/2 | -2 |
| 1/3 | -3/4 |
| 1/4 | -0.228155493653962 +1.115142508039937i |
| 1/5 | -0.154724526314600 +1.031047184160779i |
| 1/6 | i |
| 1/10 | 0.384063957 +0.666805123i |
See
* here
* Y-C Chen, T Kawahira, H-L Li, and J-M Yuan FAMILY OF INVARIANT CANTOR SETS AS ORBITS OF DIFFERENTIAL EQUATIONS. II: JULIA SETS (8475K, PDF) Dec 6, 09
The Wolf Jung test : The external parameter rays for angles (in turns)
* 1/7 (period 3)
* 321685687669320/2251799813685247 (period 51)
* 321685687669322/2251799813685247 ( period 51 )
Angles differ by about $10^{-15}$, but the landing points of the corresponding parameter rays are about 0.035 apart.

# References

* ↑ wikipedia : External ray
* ↑ gitlab repo by adam majewski: parameter external angle
* ↑ Generalizations of Douady's magic formula by Adam Epstein, Giulio Tiozzo
* ↑ fractalforums.org : generalizations-of-douadys-magic-formula
* ↑ Subdivision rule constructions on critically preperiodic quadratic matings by Mary Wilkerson
* ↑ wikipedia : Line segment
* ↑ fractalforums.org: help-with-ideas-for-fractal-art
* ↑ by Claude Heiland-Allen
* ↑ fractalforums : pathfinding-in-the-mandelbrot-set
* ↑ Family of Invariant Cantor Sets as orbits of Differential Equations II: Julia Sets by Yi-Chiuan Chen, Tomoki Kawahira, Juan-Ming Yuan
* ↑ An algorithm to draw external rays of the Mandelbrot set by Tomoki Kawahira
* ↑ fieldlines by I Quilez
* ↑ fractalforums.org: help-with-ideas-for-fractal-art
* ↑ An algorithm to draw external rays of the Mandelbrot set by Tomoki Kawahira
* ↑ Drawing external ray using Newton method ( described by T Kawahira)
* ↑ Wolf Jung's test for precision of drawing parameter rays
* Wolf Jung old page in web archive
* How to handle branch cuts in Böttcher functions
* Wolfram Research (2014), MandelbrotSetBoettcher, Wolfram Language function, https://reference.wolfram.com/language/ref/MandelbrotSetBoettcher.html.
* algorithms for computing angles in the mandelbrot set - Douady