Skip to content
Snippets Groups Projects
Select Git revision
  • master default
1 result

parameter_external_angle

  • Clone with SSH
  • Clone with HTTPS
    • adammaj1's avatar
    src
    Adam Majewski authored
    a79ee659
    History
    adammaj1's avatar a79ee659
    History
    Name Last commit Last update
    images
    src
    LICENSE
    README.md
    README.md

    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:

    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