Skip to content

adammaj1/parameter_external_angle

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 
 
 
 
 
 
 
 
 

Repository files navigation

TOC

Which points of the parameter plane have external angle ?

Points which have the external angle:

  • point of the boundary of the Mandelbrot set ( landing point of the external ray ). It can have one or more external angles : biaccesible, triaccesible , ...
  • point of the exterior of the Mandelbrot set ( only ona external angle)

Points of the interior of Mandelbrot set ( interior of Mandelbrot set hyperbolic components) do not have external angle, but have internal angle.

How to describe external angle ( proper fraction)?

Tasks related with external angle

  • compute external angle of the point c
  • draw external angle thru the point c from the exterior
  • find which external rays land on the point c from the boundary
  • draw external angle from the infinity toward boundary
  • draw/colour the parameter plane using external angle for the gradient

How to compute the external angle?

Methods for the boundary points:

Methods for the exterior points:

trace external ray outwards on the parameter plane and compute argument of last point

The lines of constant phase are exactly what is referred to as the Douady-Hubbard 'external rays'. With a tiny bit of math, its easy to see that these lines of constant phase are exactly perpendicular to the equipotential lines.

Linas Vepstas

draw ... the external rays in such a way that they are perpendicular to the escape lines

M. Romera et all in A Method to Solve the Limitations in Drawing External Rays of the Mandelbrot Set

Two methods:

  • "Choose a fixed step size, and find the direction such that a step of that size increases/decreases f the most."
  • "Choose a fixed increase in f, and find the direction such that it takes the shortest step to increase f by that amount."

Links:

series expansion formula for computing external angle

$arg_M(c) = arg(c) + \sum_{n=1}^\infty \left( \frac{1}{2^n}*arg \left( \frac{f_c^n(c)}{f_c^n(c)-c} \right ) \right ) $

The principal value of the argument is the unique angle with -π<arg(z)≤π. This definition is used because a function should be single-valued. However, a discontinuity is introduced artificially when a point z is crossing the negative real axis, say going from i through -1 to -i: its argument should go from π/2 through π to 3π/2, but the principal value goes from π/2 to π, jumps to -π, and goes to -π/2.

Wolf Jung

Note that result ot the arg function can have positive or negative sign so sum can be greater or lower then arg(c) !

whole set using palette colors

whole set using gray colors

Here one can see errors in computing : compare with Stripe Average Coloring

zoom of wake 1/3

zoom of wake 1/3 with binary decomposition

One can see on the binary decomposition image that errors are in the chaotic region, where "our image looks noisy and grainy near the boundary of the Mandelbrot set" (Claude Heiland-Allen )

Code:

/*
 Fc(z) = z*z + c
 z= x+y*i
 c= a+b*i
 This function computes external argument of point C in turns
 for mandelbrot set for Fc(z)= z*z + c
 external argument = Arg(Phi(c))
 1 [turn] = 360 [degrees] = 2* M_PI [radians]
 this function is based on function mturn from mbrot.cpp 
 from old (probably 5.2 of May 17, 2008) version of program mandel by Wolf Jung
 http://www.iram.rwth-aachen.de/~jung/indexp.html 

Already checked that escaping. Requires z = c instead of z = 0
 http://fraktal.republika.pl/cpp_argphi.html
 algorithm: http://www.mndynamics.com/indexp.html#XR

this formula will be valid not only for large |z| 

*/
double mturn(double complex c)
{ 
  int j; 
  int jMax = 100;
  
  double s = 1.0; // = 1/(2^n)
  double dr = 1.0/2.0; 
  double theta; 
  // z-c 
  double u; // creal(z-c)
  double v; // cimag(z-c)
  
  double r; //   r here means r^2 = cabs(z)* cabs(z)
  
  // z = x+y*I 
  double x; 
  double y;
  
  // c = cx + cy*I        
  double cx = creal(c);
  double cy = cimag(c);        
  
  
  
  
  // Requires z = c instead of z = 0
  x = cx;
  y = cy; 
  
  // theta = arg(z)
  theta = atan2(y, x);
  
  // compute the sum 
  for (j = 1; j < jMax; j++)
  { 
    s *= dr;
    // z=Fc(z)
    double temp = x*x - y*y + cx;
    y = 2*x*y + cy;
    x = temp; 
    
    // r here means r^2 = cabs(z)* cabs(z) 
    r = x*x  + y*y; // 
    if (r < .0001) return -7; // ?
    
    
    // z-c
    u = x - cx; 
    v = y - cy;
    
    
    
    // atan2 here is computing  arg(z/(z-c))
    // s*atan2 means : theta/(2^n)
    // theta += is summation
    theta += s * atan2(u*y - v*x, u*x + v*y);
    //
    if (r/s > 1e25) break; // prevent -nan
  }
  
  
  //if (r < 1000) return -6; // ? lazy escaping ???, 
  //
  theta *= (.5/M_PI); // convertion to turns 
  //
  theta -= floor(theta); // modulo 1 
  
  return theta;
  
  
}//mturn

atan2(y/x) = Returns the principal value of the arc tangent of y/x, expressed in radians

How to compute arg(z/(z-c)) in a fast way?

simlify z/(z-c)

z = x+y*I
arg(z) = atan2(y,x) 
z-c = u+v*I // b 
arg(z/(z-c)) = arg(z/b)

Here is Maxima CAS code:

(%i2) z:x+y*%i;
(%o2)            	%i y + x
(%i3) b:u+v*%i;
(%o3)            	%i v + u
(%i4) creal(z/b); 
	                      %i y + x
(%o4)                   creal(--------)
                              %i v + u
(%i5) realpart(z/b);
                              v y + u x
                              ---------
                               2    2
                              v  + u
(%i6) imagpart(z/b);
                             u y - v x
                             ---------
                              2    2
                             v  + u


Denote :

$z/b = rn/rd + I*in/id$

notice that

$atan2(in/id, rn/rd) = atan2(in, rn)$

Now one can skip:

  • 2 divisions : in/id and rn/rd
  • 2 multiplications : $v*v$ and u*u
  • 1 addition $v^2+u^2$

Files:

Links:

trace external ray outwards on the parameter plane

"mostly adopted a calculation method (of the external angle is) to trace the curve of the external ray"

Souichiro-Ikebe ( automatic translation)

Testing shows the original atan2() is only accurate to around 16 bits, (so) bit collection when passing dwell bands is much more accurate.

double externalAngle(...) {
...
	return (std::atan2(cy,cx));
}

This gets you the angle in only double-precision, but using double precision floating point throughout it's possible to get the external angle in much higher precision

  • the trick is to collect bits from the binary representation of the angle as you cross each dwell band
  • whether the final iterate that escaped has a positive or negative imaginary part determines if the bit is 0 or 1 respectively, see binary decomposition colouring

You need to trace a ray outwards, which means using different C values, and the bits come in reverse order, first the deepest bit from the iteration count of the start pixel, then move C outwards along the ray (perhaps using the newton's method of mandel-exray.pdf in reverse), repeat until no more bits left. you move C a fractional iteration count each time, and collect bits when crossing integer dwell boundaries

it is asymptotically too slow to be practical: $O(n^2)$ where n is the sum of the preperiod and period of the external angle.

Claude Heiland-Allen

See also:

Code:

  • tavis.cpp - compute external angle of point cx, cy

argument of the Boettcher coordinate

One can use argument of Boettcher coordinate for computing external argument (angle).

Mathematica Boettcher function

ArrayPlot[Table[Sin[100 Arg@MandelbrotSetBoettcher[x + I*y]], {y, -1, 1, .01}, {x, -2.6, .5, .01}], ImageSize->Full]

the result:

Douady and Hubbard method for c near the real axis

Douady A. (1986) Julia Sets and the Mandelbrot Set. In the book : The Beauty of Fractals. Springer, Berlin, Heidelberg, Print ISBN 978-3-642-61719-5

Douady and Hubbard found a simple method for computing external angles for values of c outside of M and near the real axis. Call such an angle 2Pi*Ray, where 0 <= Ray < 1.
The number Ray can be written as a binary decimal, i.e, as a sequence of zeroes and ones. To find it, consider the sequence {Arg[c], Arg[c^2 +c], Arg[(c^2 + c)^2 + c], ...}.

We replace Arg[z] by

  • 0 if 0 <= Arg[z] < Pi,
  • 1 otherwise

Here is some Mathematica code for this.

    c = -.75 +.0001*I; 
    z = 0;
    Do[z = z^2 + c; Print[Abs[Floor[Arg[z]/Pi]]], {n, 1, 10}]

This produces the sequence {0, 1, 0, 1, 0, 1, 0, ...} which is the binary expansion for 1/3
For c = -.75 - .0001*I produces {1, 0, 1, 0, 1, 0, 1, ...} which is the binary expansions for 2/3.
The point c0 = -.75 is the root of the period 2 bud. There are two rays leading inward to it, one coming from above and one from below. The two values of c we have chosen lie on or very near these two rays.

Douglas C. Ravenel

Files:

  • douady.c - c file wich checks Douady-Hubbard method
  • morse.mac - batch file for Maxima cas which computes upper angles of external rays which land on the roots of the period doubling cascade on the real axis

Stripe Average Coloring (or Method) = SAM or SAC

Links:

whole set

zoom of wake 1/3

Files:

See also

# escape time field line algorithm in Ultra Fractal by  Chris Thomasson
ct_test_mandelbrot_escape_field {
init:
z = #pixel
zp = z
g = 200.0

loop:
zp = z
z = z^2 + #pixel
bailout:
! (((real(z) / real(zp)) > g) && ((imag(z) / imag(zp)) > g))
default:
title = "CT Test Mandelbrot Escape Field"
}

Key words

  • discrete local complex dynamics
  • complex quadratic polynomial
  • basin of attraction of infinity
  • basin of attraction of superattracting fixed point

technical notes

GitLab uses:

git ( gitlab)

cd existing_folder
git init
git remote add origin git@gitlab.com:adammajewski/parameter_external_angle.git
git add .
git commit -m "Initial commit"
git push -u origin master

local repo : ~/c/mandel/p_e_angle

Subdirectory

mkdir images
git add *.png
git mv  *.png ./images
git commit -m "move"
git push -u origin master

then link the images:

![](./images/n.png "description") 
gitm mv -f