Jump to content

File:InfoldingSiegelDisk1over3.gif

Page contents not supported in other languages.
This is a file from the Wikimedia Commons
From Wikipedia, the free encyclopedia

Original file (999 × 999 pixels, file size: 55 KB, MIME type: image/gif, looped, 10 frames, 10 s)

Summary

Description
English: Infolding Siegel Disk near 1/3. It is a numerical aproximation so some shapes are not precise
Date
Source I have made it with significant help of Claude Heiland-Allen and Wolf Jung. This image is based on the idea taken from image by Arnaud Chéritat[1].
Author Adam majewski

Licensing

I, the copyright holder of this work, hereby publish it under the following license:
w:en:Creative Commons
attribution share alike
This file is licensed under the Creative Commons Attribution-Share Alike 4.0 International license.
You are free:
  • to share – to copy, distribute and transmit the work
  • to remix – to adapt the work
Under the following conditions:
  • attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
  • share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.

cpp source code

 
/*
 
errors : why graphic file has 999 pixel not 1000 ? 

  c console program
  It can be compiled and run under Linux, windows, Mac 
  It needs gcc
 
  --------------------------------------
  draws critical orbit for f(z)=z*z+c 


  forward 
  backward


===========================
Q: 
Let's take 
t =  0.3332223278292314 

find bifurcation point from period 1

so


c = -0.124396035791842  +0.649518736914556 i    period = 10000


then go to dynamic plane
increase iteration number to see better Julia set
then draw better critical orbit


then go back to critical point z=0

Then use
AA BAA AAA........

to stay on critical orbit ( and to draw rare points of critical orbits ) 

Is it possiible to find other sequences of inverse iterations ( A, B)
to draw rare points of such critical orbits ?
------------ Answer -----------------------
No,  there is only one backwards orbit contained in the boundary of the
Siegel disk.  Ideally it should be AAAAAAAAAAAAAAAAAAA... , but since
the commands  a  and  b  work only approximately, it might be different.

  Best regards,

  Wolf

  
 
 
 
  ------------------------------------------
  one can change :
  
  - n 
  - iSide ( width of image = iXmax = (iSide) 
  - NrOfCrOrbitPoints = ; // check rang of type for big numbers : iMax, i, NrOfCrOrbitPoints


  %lld and %llu for print long long int
  -----------------------------------------
  1.pgm file code is  based on the code of Claudio Rocchini
  http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
  create 8 bit color graphic file ,  portable gray map file = pgm 
  see http://en.wikipedia.org/wiki/Portable_pixmap
  to see the file use external application ( graphic viewer)
  I think that creating graphic can't be simpler
  ---------------------------
  2. first it creates data array which is used to store color values of pixels,
  fills tha array with data and after that writes the data from array to pgm file.
  It alows free ( non sequential) acces to "pixels"
  -------------------------------------------
  Here are 4 items : 
  1. complex plane Z with points z = zx+zy*i ( dynamic plane and world coordinate )
  2. virtual 2D array indexed by iX and iYmax ( integer or screen coordinate )
  3. memory 1D array data  indexed by i =f(iX,iY)
  4. pgm file 
 
 
 
  Adam Majewski   fraktal.republika.pl 
 
 
  to compile : 
  gcc d.c -lm -Wall -march=native
  to run ( Linux console) :
  time ./a.out



  convert -version
 
  convert data0.pgm -convolve "-1,-1,-1,-1,8,-1,-1,-1,-1"  -threshold 5% -negate data0.png
  convert data0.pgm -convolve "0,1,0,1,1,1,0,1,0" -threshold 5% -negate data0c.png
  convert data0.pgm -convolve "-0.125,0,0.125,  -0.25,0,0.25,  -0.125,0,0.125" -threshold 5% -negate data0g.png

  convert data0.pgm -edge 3 -negate data0e.png
  convert data0.pgm -edge 3 -negate data0e.png
  http://www.imagemagick.org/Usage/transform/#vision
  "As you can see, the edge is added only to areas with a color gradient that is more than 50% white! I don't know if this is a bug or intentional, but it means that the edge in the above is located almost completely in the white parts of the original mask image. This fact can be extremely important when making use of the results of the "-edge" operator. For example if you are edge detecting an image containing an black outline, the "-edge" operator will 'twin' the black lines, producing a weird result." 

  convert data0.pgm -negate -edge 3 data0n.png
  convert data0n.png -edge 3 -negate data0e.png


  http://unix.stackexchange.com/questions/299218/imagemagick-how-to-thicken-lines

  convert 4.pgm -negate 4n.pgm
  convert 4n.pgm -morphology Dilate Octagon 4nf.pgm
  convert 4n.pgm -morphology Thicken '3x1+2+0:1,0,0' 4nfb.pgm
  convert 4n.pgm -morphology Thicken ConvexHull 4nfc.pgm
  convert 4nf.pgm -negate 4thick.pgm


  --------------------
  check if curve is closed : 
  if x(1) == x(end) && y(1) == y(end)
  % It's closed
  else
  % It's not closed
  end


  ---------------------------------------------------------------------------------

  http://www.scholarpedia.org/article/File:InfoldingSiegelDisk.gif
  1,000 × 1,000 pixels, file size: 91 KB, MIME type: image/gif, looped, 9 frames, 12s)
  The Siegel disks have been translated in the plane so that the critical point is always at the same point on the screen (near the top). 


  for n :
  0,1   I have used  1.0e5 = pow(10,5) points 
  n= 2                    1.0e6 = pow(10,6)
  n = 3                   1.0e7 = pow(10,7)
  n = 4                   1.0e8 = pow(10,8) // good result
  n= 5                    1.0e9 = pow(10,9) // not good 
 
  For n=5 I have to try pow(10,12).
 
  Do you use such high number of iterations or other method ?

  I think in this particular case, I iterated a lot. However, in some other pictures, I used an acceleration method, an approximation of a high iterate of f.



  










*/

#define BOUNDS_CHECKS
// uncomment next line for tiny speedup but less safety (possible invalid memory writes)
//#undef BOUNDS_CHECKS

#ifdef __cplusplus

#include <complex>

#ifdef USE_QD_REAL
#define QD_INLINE
#include <qd/qd_real.h>
typedef qd_real Real;
#define get_double(z) ((z).x[0])
#define pi (Real("3.141592653589793238462643383279502884197169399375105820974944592307816406286208"))
#else
#ifdef USE_DD_REAL
#define QD_INLINE
#include <qd/dd_real.h>
typedef dd_real Real;
#define get_double(z) ((z).x[0])
#define pi (Real("3.141592653589793238462643383279502884197169399375105820974944592307816406286208"))
#else
typedef double Real;
#define get_double(z) (z)
#define pi (3.141592653589793238462643383279502884197169399375105820974944592307816406286208)
static inline Real sqr(Real x) { return x * x; }
static inline Real mul_pwr2(Real x, double y) { return x * y; }
#endif
#endif

typedef std::complex<Real> Cplx;

#define creal(c) (std::real(c))
#define cimag(c) (std::imag(c))
#define I (Cplx(0,1))

#else

#include <complex.h>

typedef double Real;
typedef double _Complex Cplx;

#endif


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

int n; 
#define NMAX 13
#define L  (NMAX +1)

#define unlikely(x) (__builtin_expect(!!(x),0))
 
/* iXmax/iYmax =  */
const int iSide = 1000;
int iXmax ; /* width of image in pixels = (15*iSide); */
int iYmax ;
int iLength ;

/* dynamic 1D arrays for colors ( shades of gray )  */
  
unsigned char *data;
unsigned char *data3;


/* */
double ZxMin = -0.75;
double ZxMax = 0.35;
double ZyMin  = -0.15;
double ZyMax = 0.85;
/* (ZxMax-ZxMin)/(ZyMax-ZyMin)= iXmax/iYmax  */
 
 
double PixelWidth ;
double PixelHeight ;
double invPixelWidth ;
double invPixelHeight ;
 
unsigned int period=1;
unsigned int iPeriodChild=3;

// unsigned int m;
// check rang of type for big numbers : iMax, i, NrOfCrOrbitPoints
unsigned long long int NrOfCrOrbitPoints ;
 
double InnerRadius= 0.05997604;
 
Real radius = 1.0;
Real tt(int n)
{
  // t(n) = [0;3,  10^n, golden_ratio]
  Real phi = (sqrt(Real(5)) + 1) / 2;
  Real tenn = pow(Real(10), n);
  return 0 +1/( 3  +1/( tenn +1/( phi )));
}

/* fc(z) = z*z + c */
/*   */
Cplx C; 
Real Cx; /* C = Cx + Cy*i   */
Real Cy ;
 
Cplx alfa; 
Real Ax;
Real Ay;

Cplx z3a, z3b, z3c;
 
/* colors */
const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
const int iExterior = 150; /* exterior of Julia set */
const int iJulia = 0; /* border , boundary*/
const int iJuliaInv = 50; /* border , boundary*/
const int iInterior = 200;
const int iBackground = 255; 
 
 
/* ----------------------- functions ---------------------------------------- */


/* find c in component of Mandelbrot set 
 
   uses code by Wolf Jung from program Mandel
   see function mndlbrot::bifurcate from mandelbrot.cpp
   http://www.mndynamics.com/indexp.html

*/
Cplx GiveC(Real InternalAngleInTurns, Real InternalRadius, unsigned int Period)
{
  //0 <= InternalRay<= 1
  //0 <= InternalAngleInTurns <=1
  Real t = InternalAngleInTurns *2*pi; // from turns to radians
  Real R2 = InternalRadius * InternalRadius;
  //Real Cx, Cy; /* C = Cx+Cy*i */
  switch ( Period ) // of component 
    {
    case 1: // main cardioid
      Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4; 
      Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4; 
      break;
    case 2: // only one component 
      Cx = InternalRadius * 0.25*cos(t) - 1.0;
      Cy = InternalRadius * 0.25*sin(t); 
      break;
      // for each iPeriodChild  there are 2^(iPeriodChild-1) roots. 
    default: // higher periods : to do, use newton method 
      Cx = 0.0;
      Cy = 0.0; 
      break; }

  return Cx + Cy*I;
}


/*

  http://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings
  z^2 + c = z
  z^2 - z + c = 0
  ax^2 +bx + c =0 // ge3neral for  of quadratic equation
  so :
  a=1
  b =-1
  c = c
  so :

  The discriminant is the  d=b^2- 4ac 

  d=1-4c = dx+dy*i
  r(d)=sqrt(dx^2 + dy^2)
  sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx +- sy*i

  x1=(1+sqrt(d))/2 = beta = (1+sx+sy*i)/2

  x2=(1-sqrt(d))/2 = alfa = (1-sx -sy*i)/2

  alfa : attracting when c is in main cardioid of Mandelbrot set, then it is in interior of Filled-in Julia set, 
  it means belongs to Fatou set ( strictly to basin of attraction of finite fixed point )

*/
// uses global variables : 
//  ax, ay (output = alfa(c)) 
Cplx GiveAlfaFixedPoint(Cplx c)
{
  Real dx, dy; //The discriminant is the  d=b^2- 4ac = dx+dy*i
  Real r; // r(d)=sqrt(dx^2 + dy^2)
  Real sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i
  Real ax, ay;
 
  // d=1-4c = dx+dy*i
  dx = 1 - 4*creal(c);
  dy = -4 * cimag(c);
  // r(d)=sqrt(dx^2 + dy^2)
  r = sqrt(dx*dx + dy*dy);
  //sqrt(d) = s =sx +sy*i
  sx = sqrt((r+dx)/2);
  sy = sqrt((r-dx)/2);
  // alfa = ax +ay*i = (1-sqrt(d))/2 = (1-sx + sy*i)/2
  ax = 0.5 - sx/2.0;
  ay =  sy/2.0;
 

  return ax+ay*I;
}
 

 
 
inline unsigned int f(unsigned int iX, unsigned int iY)
/* 
   gives position of point (iX,iY) in 1D array  ; uses also global variables 
   it does not check if index is good  so memory error is possible 
*/
{return (iX + iY*iXmax );}
 
 
 
inline int DrawPoint( Real Zx, Real Zy, unsigned char data[], int iColor)
{
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int index; /* index of 1D array  */ 
 
#ifdef BOUNDS_CHECKS
  //if (unlikely(Zx < ZxMin || ZxMax  < Zx || Zy < ZyMin || ZyMax < Zy)) {  printf("   point z = %f , %f out of bounds  \n", get_double(Zx), get_double(Zy)); return -1; } 
#endif
  
  iX = (int)get_double((Zx-ZxMin)*invPixelWidth);
  iY = (int)get_double((ZyMax-Zy)*invPixelHeight); // reverse Y axis
  index = iX + iY*iXmax;//f(iX,iY);
  
  
  data[index] = iColor;  /* draw */
  return 0;
}
 


int TestCriticalOrbit(unsigned long long int iMax )
{
  unsigned long long int i; /* nr of point of critical orbit */
  Real Zx,Zy, tmp;
   
 
  /* critical point z = 0 */
  Zx = 0.0;
  Zy = 0.0;
  //
  ZxMin = 0.0;
  ZxMax = 0.0;
  ZyMin = 0.0;
  ZyMax = 0.0;
  
  /* forward orbit of critical point  */
  for (i=1;i<=iMax ;i++)
    {
      /* f(z) = z*z+c */
      tmp = mul_pwr2(Zx*Zy, 2) + Cy;
      Zx = sqr(Zx) - sqr(Zy) + Cx;
      Zy = tmp;
       
 
      
      // if (Zx2+Zy2>4) { printf("   point z = %f , %f escapes \n",Zx, Zy); break;}
      if (Zx>ZxMax) ZxMax=get_double(Zx);
      if (Zx<ZxMin) ZxMin=get_double(Zx);
      if (Zy>ZyMax) ZyMax=get_double(Zy); 
      if (Zy<ZyMin) ZyMin=get_double(Zy);         

    }

  printf(" ZxMin = %.16f  ZxMax = %.16f \n", ZxMin, ZxMax);       
  printf(" ZyMin = %.16f  ZyMax = %.16f \n", ZyMin, ZyMax);  
  return 0;
 
}


/* mndyncxmics::root from mndyncxmo.cpp  by Wolf Jung (C) 2007-2014. */

// input = x,y
// output = u+v*I = sqrt(x+y*i) 
Cplx
GiveRoot (Cplx z)
{
  Real x = creal (z);
  Real y = cimag (z);
  Real u, v;

  v = sqrt (x * x + y * y);

  if (x > 0.0)
    {
      u = sqrt (0.5 * (v + x));
      v = 0.5 * y / u;
      return u + v * I;
    }
  if (x < 0.0)
    {
      v = sqrt(Real(0.5 * (v - x)));
      if (y < 0.0)
	v = -v;
      u = 0.5 * y / v;
      return u + v * I;
    }
  if (y >= 0.0)
    {
      u = sqrt(Real(0.5 * y));
      v = u;
      return u + v * I;
    }

  u = sqrt (Real (-0.5 * y));
  v = -u;
  return u + v * I;
}

// from mndlbrot.cpp  by Wolf Jung (C) 2007-2014. part of Madel 5.12 
// input : c, z , mode
// c = cx+cy*i where cx and cy are global variables defined in mndynamo.h
// z = x+y*i
// 
// output : z = x+y*i
Cplx InverseIteration(Cplx z, Cplx c)
{
  Real x = creal (z);
  Real y = cimag (z);
  Real cx = creal (c);
  Real cy = cimag (c);

  // f^{-1}(z) = inverse with principal value
  if (cx * cx + cy * cy < 1e-20)
    {
      z = GiveRoot (x - cx + (y - cy) * I);	// 2-nd inverse function = key b 
      //if (mode & 1) { x = -x; y = -y; } // 1-st inverse function = key a   
      return -z;
    }

  //f^{-1}(z) =  inverse with argument adjusted
  Real u, v;
  Cplx uv;
  Real w = cx * cx + cy * cy;

  uv = GiveRoot (-cx / w - (cy / w) * I);
  u = creal (uv);
  v = cimag (uv);
  //
  z = GiveRoot (w - cx * x - cy * y + (cy * x - cx * y) * I);
  x = creal (z);
  y = cimag (z);
  //
  w = u * x - v * y;
  y = u * y + v * x;
  x = w;

  //if (mode & 1) // mode = -1
  //  { x = -x; y = -y; } // 1-st inverse function = key a

  return x + y * I;		// 2-nd inverse function = key b 
}


 int DrawInverseOrbit(unsigned long long int iMax,  unsigned char A[], unsigned char iColor ){

    Cplx Z;

/* critical point z = 0 */
 
 
 Z = -  InverseIteration(Cplx(0.0,0.0), C); // A
 DrawPoint(get_double(creal(Z)),get_double(cimag(Z)),A, iColor);
 
 Z = -  InverseIteration(Z, C); // A
 DrawPoint(get_double(creal(Z)),get_double(cimag(Z)),A, iColor);

 Z =   InverseIteration(Z, C); // B
 DrawPoint(get_double(creal(Z)),get_double(cimag(Z)),A, iColor);
 for (unsigned long long int i = 0; i < iMax; ++i){
 Z = -  InverseIteration(Z, C); // A
 DrawPoint(get_double(creal(Z)),get_double(cimag(Z)),A, iColor);
 }
 


return 0; 

}


// forward iteration
// Cx and Cy = global
Cplx Forward(Cplx z){
  Real Zx,Zy, tmp;
  Zx = creal(z);
  Zy = cimag(z);

	  /* f(z) = z*z+c */
	  tmp = mul_pwr2(Zx*Zy, 2) + Cy;
	  Zx = sqr(Zx) - sqr(Zy) + Cx;
	  Zy = tmp;
   return Cplx(Zx,Zy);

}

int DrawForwardOrbit(Cplx Z, unsigned long long int iMax,  unsigned char A[], unsigned char iColor  )
{
  //unsigned long long int iMax10000 = iMax / 10000;
  unsigned long long int i; /* nr of point of critical orbit */
  Real Zx,Zy, tmp;
  int IsGood;  
 
  /* Z0 */
  Zx = creal(Z);
  Zy = cimag(Z);
  DrawPoint(Zx,Zy,A,iColor);
  
   
      for (i=1;i<=iMax ;i++)
	{
	  /* f(z) = z*z+c */
	  tmp = mul_pwr2(Zx*Zy, 2) + Cy;
	  Zx = sqr(Zx) - sqr(Zy) + Cx;
	  Zy = tmp;
       
 
      
	  // if (Zx2+Zy2>4) { printf("   point z = %f , %f escapes \n",Zx, Zy); break;}
      
	 // IsGood =
          DrawPoint(Zx,Zy,A,iColor); /* draws critical orbit */
#ifdef BOUNDS_CHECKS
	//  if (unlikely(IsGood<0)) { printf ("i = %llu \n", i); return 1;}
#endif
	}
    
    
   
  return 0;
 
}

 
/*

  check rang of type for big numbers : iMax, i, NrOfCrOrbitPoints
  x := x^2 - y^2 + cx
  y := 2 x y + cy
*/


int DrawCriticalOrbit(unsigned long long int iMaxInverse, unsigned long long int iMaxForward, unsigned char A[] )
{
  
  unsigned long long int i; /* nr of point of critical orbit */
  Cplx Z0;
  Cplx Z; 
 
   
  /* critical point z = 0 */
  Z0 = Cplx(0.0,0.0);
  //printf("   point z = %f , %f   \n", get_double(creal(Z)), get_double(cimag(Z)));
  printf("%dC",n); // symbolic sequence starting from n , then  crritical ,  

  DrawForwardOrbit(Z0, iMaxForward, A, iJulia);



  // first
  Z =  InverseIteration(Z, C); 
  if (cimag(Z) >0.0) // manual adjust 
                  {printf("B"); }
             else {Z = -Z; printf("A"); } 
  //printf("   point z = %f , %f   \n", get_double(creal(Z)), get_double(cimag(Z)));
  DrawPoint( get_double(creal(Z)), get_double(cimag(Z)),A,iJuliaInv);


  // second
  //DrawForwardOrbit(Z, 2, A);
  Z =  InverseIteration(Z, C); 
       if ( get_double(cimag(Z)) >0.0) 
                  {printf("B"); }
             else {Z = -Z; printf("A"); } 
  //printf("   point z = %f , %f   \n", get_double(creal(Z)), get_double(cimag(Z)));
 // DrawForwardOrbit(Z, 2, A);
  DrawPoint( get_double(creal(Z)), get_double(cimag(Z)),A,iJuliaInv);




  // third 
  Z =  InverseIteration(Z, C); 
       if ( get_double(cimag(Z)) >0.0) 
                  {printf("B"); }
             else {Z = -Z; printf("A"); } 

  //printf("   point z = %f , %f   \n", get_double(creal(Z)), get_double(cimag(Z)));



   for (i=1;i<=iMaxInverse ;i++){
     //DrawForwardOrbit(Z,2, A);
       DrawPoint( get_double(creal(Z)), get_double(cimag(Z)),A,iJuliaInv);
       Z =  InverseIteration(Z, C); 
       if ( get_double(cimag(Z)) >0.0) 
                  {printf("B"); }
             else {Z = -Z; printf("A"); } 
      // if ((get_double(creal(Z)) < ZxMin || ZxMax  < get_double(creal(Z)) || get_double(cimag(Z)) < ZyMin || ZyMax < get_double(cimag(Z)))) {  printf("   point z = %f , %f out of bounds  \n", get_double(creal(Z)), get_double(cimag(Z))); return -1; } 
    }
   printf("\n"); // symbolic sequence 
  return 0;
 
}



/*
The distance between (x1, y1) and (x2, y2) is given by:
https://en.wikipedia.org/wiki/Distance

*/

Real GiveD(Real x1, Real y1, Real x2, Real y2 ){

return sqrt((x1-x2)*(x1-x2) +(y1-y2)*(y1-y2));


}

int MarkSets(unsigned long long int iterMax,  unsigned char A[] )
{
  int ix, iy; // pixel coordinate 
  Real Zx, Zy; //  Z= Zx+ZY*i;
  Real Zx2, Zy2;
  Real d; 
  
  int i; /* index of 1D array */
  unsigned int iter;
  for(iy=0;iy<=iYmax;++iy) 
    {
      //Zy =ZyMax - iy*PixelHeight;
      for(ix=0;ix<=iXmax;++ix) 
	{
          i = f(ix, iy); /* compute index of 1D array from indices of 2D array */
          if ( A[i] != iJulia) {

	  // from screen to world coordinate
          Zy =ZyMax - iy*PixelHeight; 
	  Zx = ZxMin + ix*PixelWidth ;
	  
          Zx2=Zx*Zx;
          Zy2=Zy*Zy;

          for (iter=1;iter<=iterMax ;iter++)
	{
	  /* f(z) = z*z+c */
	  Zy=2*Zx*Zy + Cy;
          Zx=Zx2-Zy2 +Cx;
          Zx2=Zx*Zx;
          Zy2=Zy*Zy;
       
        if (Zx2+Zy2>4) {A[i]=iExterior; break;}
        d = GiveD(Zx, Zy, Ax, Ay);
        if ( d< InnerRadius) {A[i]=iInterior; break;}
        

       } //iter
        
      } // if ( A[i] != iJulia)
	  

	} // ix
    } //iy
   
    
   
  return 0;
 
}





/* 
   close the curve by filing gaps by streight lines
   curve is the simple closed 2d plane curve 

*/
int CloseTheCurve( unsigned char A[] ){
  (void) A;
  return 0;
}


int ClearArray( unsigned char A[] )
{
  int index; /* index of 1D array  */
  for(index=0;index<iLength-1;++index) 
    A[index]=iBackground;
  return 0;
}


int ClearArray3( unsigned char A[] )
{
  int index; /* index of 1D array  */

  int k ; //
  for(index=0;index<iLength-1;++index) {

    k = index * 3;
    //printf("index = %d k = %d \n", index,k);
    A[k]   = iBackground;
    A[k+1] = iBackground;
    A[k+2] = iBackground;
   }
  return 0;
}


// https://commons.wikimedia.org/wiki/File:Eight_point_neighborhood.svg
int GiveNeighbours(unsigned char S[], int i){

int neighbours =0;
int iX, iY;
div_t result; // ang. result = iloraz 
  
 // i = iX + iY*iXmax;
 // compute iX and iY from i
 
 result = div(i, iXmax);
 iY = result.quot;
 iX = result.rem;

 if ( iX<1 || iY <1 || iX>iXmax-1 || iY> iYmax-1 ) {printf("out of bounds in GiveNeighbours\n"); return -1;}
 
 // iY-1
  if (S[f(iX-1,iY-1)]==iJulia) neighbours+=1;
  if (S[f(iX,iY-1)]==iJulia) neighbours+=1;
  if (S[f(iX+1,iY-1)]==iJulia) neighbours+=1;
 // iY
  if (S[f(iX-1,iY)]==iJulia) neighbours+=1;
  if (S[f(iX+1,iY)]==iJulia) neighbours+=1;
 // iY+1
  if (S[f(iX-1,iY+1)]==iJulia) neighbours+=1;
  if (S[f(iX,iY+1)]==iJulia) neighbours+=1;
  if (S[f(iX+1,iY+1)]==iJulia) neighbours+=1;




return neighbours;


}



/*
 S = source = 1 bit color
 D = destination = 3 bit color
*/

int CopyArray( unsigned char S[], unsigned char D[] )
{
  int i; /* index of 1D array  */
  //int Neighbours;
  int k ; //
  for(i=0;i<iLength-1;++i) {

     k = i*3;
    if (S[i]==iJulia){

   // Neighbours = GiveNeighbours(S,i);
    
   
    //printf("index = %d k = %d \n", index,k);

    //if (Neighbours<2) // disconnected 
       {D[k]=255;
        D[k+1]=0;
        D[k+2]=0;}
    //  else { // connected
    //    D[k]=0;
    //    D[k+1]=0;
   //     D[k+2]=255;}
    } // julia


     if (S[i]==iJuliaInv){

    //Neighbours = GiveNeighbours(S,i);
    
   
    //printf("index = %d k = %d \n", index,k);

    //if (Neighbours<2) // disconnected 
       {D[k]=0;
        D[k+1]=0;
        D[k+2]=255;}
    //  else { // connected
     //   D[k]=0;
     //   D[k+1]=255;
      //  D[k+2]=255;}
    } // julia



    
    if (S[i]==iInterior){
         D[k]=iInterior;
        D[k+1]=iInterior;
        D[k+2]=iInterior;}// interior
    


    if (S[i]==iExterior){
         D[k]=iExterior;
        D[k+1]=iExterior;
        D[k+2]=iExterior;}// interior
     //
} //
   
  return 0;
}

/*
 S = source = 1 bit color
 D = destination = 3 bit color
*/

int CopyTestArray( unsigned char S[], unsigned char D[] )
{
  int i; /* index of 1D array  */
  int Neighbours;
  int k ; //
  for(i=0;i<iLength-1;++i) {

     k = i*3;
    if (S[i]==iJulia){

    Neighbours = GiveNeighbours(S,i);
    
   
    //printf("index = %d k = %d \n", index,k);

    if (Neighbours<2) // disconnected 
       {D[k]=255;
        D[k+1]=0;
        D[k+2]=0;}
      else { // connected
        D[k]=0;
        D[k+1]=0;
        D[k+2]=255;}
    } // julia


     if (S[i]==iJuliaInv){

    Neighbours = GiveNeighbours(S,i);
    
   
    //printf("index = %d k = %d \n", index,k);

    if (Neighbours<2) // disconnected 
       {D[k]=255;
        D[k+1]=0;
        D[k+2]=255;}
      else { // connected
        D[k]=0;
        D[k+1]=255;
        D[k+2]=255;}
    } // julia



    
    if (S[i]==iInterior){
         D[k]=iInterior;
        D[k+1]=iInterior;
        D[k+2]=iInterior;}// interior
    


    if (S[i]==iExterior){
         D[k]=iExterior;
        D[k+1]=iExterior;
        D[k+2]=iExterior;}// interior
     //
} //
   
  return 0;
}


// Check Orientation of image : first quadrant in upper right position
// uses global var :  ...
int CheckOrientation(unsigned char A[] )
{
  int ix, iy; // pixel coordinate 
  Real Zx, Zy; //  Z= Zx+ZY*i;
  int i; /* index of 1D array */
  for(iy=0;iy<=iYmax;++iy) 
    {
      Zy =ZyMax - iy*PixelHeight;
      for(ix=0;ix<=iXmax;++ix) 
	{

	  // from screen to world coordinate 
	  Zx = ZxMin + ix*PixelWidth ;
	  i = f(ix, iy); /* compute index of 1D array from indices of 2D array */
	  if (Zx>0 && Zy>0) A[i] = 255-A[i];   // check the orientation of Z-plane by marking first quadrant */

	}
    }
   
  return 0;
}

/*
http://rosettacode.org/wiki/Bitmap/Bresenham%27s_line_algorithm
Instead of swaps in the initialisation use error calculation for both directions x and y simultaneously:
*/
void iDrawLine( int x0, int y0, int x1, int y1, unsigned char iColor, unsigned char A[]) 
{
  int x=x0; int y=y0;
  int dx = abs(x1-x0), sx = x0<x1 ? 1 : -1;
  int dy = abs(y1-y0), sy = y0<y1 ? 1 : -1; 
  int err = (dx>dy ? dx : -dy)/2, e2;
  int index; 
  e2 = err;
  for(;;){
    //iDrawPoint(A, x, y,color);

   index = x + y*iXmax;//f(iX,iY);
   A[index] = iColor;  /* draw */




   if (x==x1 && y==y1) break;
   // if (x<0 || x>iXmax || y<0 || y > iYmax ) break; // cliping 
    
    if (e2 >-dx) { err -= dy; x += sx; }
    if (e2 < dy) { err += dx; y += sy; }
    e2 = err;
    
  }
}
/*

paritition a dynamic plane into k=PeriodChild sectors
with center = alf afixed point 
when point z from Siegel disc if iterated 
it moves from one sector to other
*/


int dDrawLine(double Zx0, double Zy0, double Zx1, double Zy1, unsigned char iColor, unsigned char A[]) 
{

  int ix0, iy0; // screen coordinate = indices of virtual 2D array 
  int ix1, iy1; // screen coordinate = indices of virtual 2D array

   // first step of clipping
   if (  Zx0 < ZxMax &&  Zx0 > ZxMin && Zy0 > ZyMin && Zy0 <ZyMax 
      && Zx1 < ZxMax &&  Zx1 > ZxMin && Zy1 > ZyMin && Zy1 <ZyMax )
   ix0= (Zx0- ZxMin)/PixelWidth; 
   iy0 = (ZyMax - Zy0)/PixelHeight; // inverse Y axis 
   ix1= (Zx1- ZxMin)/PixelWidth; 
   iy1= (ZyMax - Zy1)/PixelHeight; // inverse Y axis 
   // second step of clipping
   if (ix0 >=0 && ix0<=iXmax && ix1 >=0 && ix1<=iXmax && iy0 >=0 && iy0<=iYmax && iy1 >=0 && iy1<=iYmax )

   //
   iDrawLine(ix0,iy0,ix1,iy1, iColor, A) ;

return 0;
}

/* 
from mndlbrot.cpp  part of Mandel 5.12 by Wolf Jung (C) 2007-2014.  
 the GNU General Public License as published by the Free Software Foundation; either
   version 3, or (at your option) any later version.

int find(int sg, unsigned int preper, unsigned int per, double cx, double cy, double &x, double &y) 
{  double a = cx, b = cy, fx, fy, px, py, w; 
   unsigned int i, j;
   for (i = 0; i < 30; i++)
   {  if (sg > 0) { a = x; b = y; }
      if (!preper)
      {  if (sg > 0) { fx = 0; fy = 0; px = 0; py = 0; }
         else { fx = -x; fy = -y; px = -1; py = 0; }
      }
      else
      {  fx = x; fy = y; px = 1.0; py = 0;
         for (j = 1; j < preper; j++)
         {  if (px*px + py*py > 1e100) return 1;
            w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
            px = w; if (sg > 0) px++;
            w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
         }
      }
      double Fx = fx, Fy = fy, Px = px, Py = py;
      for (j = 0; j < per; j++)
      {  if (px*px + py*py > 1e100) return 2;
         w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
         px = w; if (sg > 0) px++;
         w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
      }
      fx += Fx; fy += Fy; px += Px; py += Py;
      w = px*px + py*py; if (w < 1e-100) return -1;
      x -= (fx*px + fy*py)/w; y += (fx*py - fy*px)/w;
   }
   return 0;
}



*/













Cplx  find( int preper, int period, Real a, Real b) 
{  
   Real x=0.0;
   Real y=0.0;
   Real  fx=0.0;
   Real  fy=0.0;
   Real px=0.0;
   Real py=0.0;
   Real w=0.0; 
   int i, j;

   Real Px = 0.0;
   Real Py = 0.0;
   Real Fx=0.0;
   Real Fy=0.0;


    for (i = 0; i < 30; i++)
   { 
      if (!preper)
      { { fx = -x; fy = -y; px = -1.0; py = 0.0; }
      }
      else
      {  fx = x; fy = y; px = 1.0; py = 0.0;
         for (j = 1; j < preper; j++)
         {  if (px*px + py*py > 1e100) return Cplx(0.0,0.0);
            w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
            px = w; 
            w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
         }
      }
      Real Fx = fx, Fy = fy, Px = px, Py = py;
      for (j = 0; j < period; j++)
      {  if (px*px + py*py > 1e100) Cplx(1.0,0.0);
         w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
         px = w;  
         w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
      }
      fx += Fx; fy += Fy; px += Px; py += Py;
      w = px*px + py*py; if (w < 1e-100) return Cplx(-1,0.0);
      x -= (fx*px + fy*py)/w; y += (fx*py - fy*px)/w;
   }
   return Cplx(x,y);
}


int DrawStar(unsigned char A[] ){

  dDrawLine(get_double(Ax), get_double(Ay),get_double(creal(z3a)), get_double(cimag(z3a)), 0, A);
  dDrawLine(get_double(Ax), get_double(Ay),get_double(creal(z3b)), get_double(cimag(z3b)), 0, A);
  dDrawLine(get_double(Ax), get_double(Ay),get_double(creal(z3c)), get_double(cimag(z3c)), 0, A);

return 0; 
}



int SaveArray2pgm(unsigned char A[], unsigned int n)
{

  FILE * fp;
  char name [20]; /* name of file */
  sprintf(name,"%u",n); /*  */
  char *filename =strcat(name,".pgm");
  const char *comment="# C= ";/* comment should start with # */
  /* save image to the pgm file  */      
  fp= fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"P5\n %s\n %u %u\n %u\n",comment,iXmax,iYmax,MaxColorComponentValue);  /*write header to the file*/
  fwrite(A,iLength,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);
  return 0;
}
 




int SaveArray2ppm(unsigned char A[], unsigned int n)
{

  FILE * fp;
  char name [20]; /* name of file */
  sprintf(name,"%u",n); /*  */
  char *filename =strcat(name,".ppm");
  //const char *comment="# C= ";/* comment should start with # */
  /* save image to the pgm file  */      
  fp= fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"P6\n %u %u\n %u\n",iXmax,iYmax,MaxColorComponentValue);  /*write header to the file*/
  fwrite(A,iLength*3,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);
  return 0;
}
 

unsigned long long int GiveNrOfCrOrbitPoints( int n){

  

  //unsigned long long int i ;

  //i = 10000000000;
  // i = fabs(get_double(7 / (2 - 7 * tt(n))));

  unsigned long long int i ; 
  switch( n )
  { case  0 : i =      100000; break;
    case  1 : i =      100000; break;  
    case  2 : i =      1000000; break;
    case  3 : i =      10000000; break;
    case  4 : i =      100000000; break;
    case  5 : i =      1000000000; break;
    case  6 : i =      10000000000; break;
    case  7 : i =      10000000000; break;
    case  8 : i =      100000000000; break;
    case  9 : i =      100000000000; break;
    case 10 : i =      1000000; break;
    case 11 : i =      1000000; break;
    default: i = 100000; break;
    }
  return i;


  return i;
}

int setup(int n){

  

  /* unsigned int iX,iY;  indices of 2D virtual array (image) = integer coordinate */
  iXmax = iSide-1; /* height of image in pixels */
  iYmax = iSide-1;
  iLength = (iSide*iSide);
 
  //
  PixelWidth =  ((ZxMax-ZxMin)/iXmax);
  PixelHeight = ((ZyMax-ZyMin)/iYmax);
  invPixelWidth = 1 / PixelWidth;
  invPixelHeight = 1 / PixelHeight;
 
 
  // 1 byte color 
  data = (unsigned char *) malloc( iLength * sizeof(unsigned char) );
  if (data == NULL )
    {
      fprintf(stderr," Could not allocate memory");
      return 1;
    }
  
          
  

  // 3 byte color
  data3 = (unsigned char *) malloc( 3*iLength * sizeof(unsigned char) );
  if (data3 == NULL )
    {
      fprintf(stderr," Could not allocate memory");
      return 1;
    }
  ClearArray( data );
  ClearArray3( data3 );

   
  //
  C = GiveC(tt(n), radius, period);
  Cx = creal(C); 
  Cy = cimag(C); 
  //
  alfa = GiveAlfaFixedPoint(C);
  Ax = creal (alfa);
  Ay = cimag (alfa);
  // 
  NrOfCrOrbitPoints = GiveNrOfCrOrbitPoints(n); //pow(10,5+n);
  //

  z3a = find(0,iPeriodChild,Cx, Cy); 
  z3b = Forward(z3a);
  z3c = Forward(z3b);

  return 0;
}


int info(int n){

  printf("   period = %d  \n", period);
  printf("   n = %d  \n", n);
  printf("   t = %.16f \n", get_double(tt(n)));
  printf("   c = %.16f , %.16f  \n",get_double(Cx), get_double(Cy));
  printf("   alfa = %.16f , %.16f  \n",get_double(Ax), get_double(Ay));
  printf(" periodic orbit : \n");
  printf("   z3a = %.16f , %.16f  \n",get_double(creal(z3a)), get_double(cimag(z3a)));
  printf("   z3a = %.16f , %.16f  \n",get_double(creal(z3b)), get_double(cimag(z3b)));
  printf("   z3a = %.16f , %.16f  \n",get_double(creal(z3c)), get_double(cimag(z3c)));

  printf("   NrOfCrOrbitPoints = %.llu  \n",NrOfCrOrbitPoints);
  printf("   Inner Radius = %.16f \n", InnerRadius);
  //
  printf("   iSide = %d  \n", iSide);
  printf("   iXmax = %d  \n", iXmax);
  printf("   iLength = %d  \n", iLength);


  return 0;
}
 






/* --------------------------------------------------------------------------------------------------------- */
 
int main(){
  
    
 
  for (n=9; n<=9; ++n){ 

    setup(n);
    //TestCriticalOrbit(NrOfCrOrbitPoints);
    /* ------------------ draw ----------------------- */
    printf(" for n = %d draw %llu points of critical orbit to array \n",n, NrOfCrOrbitPoints);       
    DrawCriticalOrbit(50, NrOfCrOrbitPoints, data);
    SaveArray2pgm(data, n);
    //
    CopyArray(data,data3);
    SaveArray2ppm(data3,n);
    //
   // MarkSets(300000, data);
   // CopyArray(data,data3);
    //SaveArray2pgm(data, 100+n);
    //SaveArray2ppm(data3,100+n);
   // printf(" save  data array to the pgm file \n");
   

    


   // SaveArray2ppm(data3,100+n);
    //CheckOrientation(data);
    //SaveArray2pgm(data, n+1000);

    //DrawStar(data);
    //SaveArray2pgm(data, 100+n);


     info(n); 

  }
  

  //
  free(data);
  free(data3);
 
 
 
 
 
  return 0;
}

Makefile

all: d_d d_dd d_qd

d_d: d.c
	gcc -o d_d  -std=c++17 -Wall -pedantic -Wextra -xc++ -O3 -march=native d.c -lm -lqd -lstdc++

d_dd: d.c
	gcc -o d_dd -std=c++17 -Wall -pedantic -Wextra -xc++ -O3 -march=native d.c -lm -lqd -lstdc++ -DUSE_DD_REAL

d_qd: d.c
	gcc -o d_qd -std=c++17 -Wall -pedantic -Wextra -xc++ -O3 -march=native d.c -lm -lqd -lstdc++ -DUSE_QD_REAL

bash source code

#!/bin/bash
 
# script file for BASH 
# which bash
# save this file as g.sh
# chmod +x g.sh
# ./g.sh
 
# for all pgm files in this directory
for file in *.pgm ; do
  # b is name of file without extension
  b=$(basename $file .pgm)
  # convert  using ImageMagic
  convert $file -negate ${b}.pgm
  echo $file
done

# for all pgm files in this directory
for file in *.pgm ; do
  # b is name of file without extension
  b=$(basename $file .pgm)
  # convert  using ImageMagic
  convert $file -morphology Thicken ConvexHull ${b}.gif
  echo $file
done

 
# convert gif files to animated gif
convert -delay 100   -loop 0 %d.gif[0-11] a11.gif
 
echo OK
# end

postprocessing

I had to plot a few extra pixels with XPaint ( Fat Bits editor) to close the curve ( there was gaps at the peripheric parts of the image). I think that images should be simple 2d closed curves Beacause of that shape of curve in the peripheral parts is not good ( seee FAQ)


In last version of image inverse iteration is used without manual plotting.

text output

 ./d_qd
   period = 1  
   n = 0  
   t = 0.2763932022500210 
   c = 0.1538380639536641 , 0.5745454151066985  
   NrOfCrOrbitPoints = 1000000  
   period = 1  
   n = 1  
   t = 0.3231874668087892 
   c = -0.0703924965263780 , 0.6469145331346999  
   NrOfCrOrbitPoints = 1000000  
   period = 1  
   n = 2  
   t = 0.3322326933513446 
   c = -0.1190170769366243 , 0.6494880316361160  
   NrOfCrOrbitPoints = 1000000  
   period = 1  
   n = 3  
   t = 0.3332223278292314 
   c = -0.1243960357918422 , 0.6495187369145560  
   NrOfCrOrbitPoints = 1000000  
   period = 1  
   n = 4  
   t = 0.3333222232791965 
   c = -0.1249395463818515 , 0.6495190496732967  
   NrOfCrOrbitPoints = 1000000  
   period = 1  
   n = 5  
   t = 0.3333322222327929 
   c = -0.1249939540657307 , 0.6495190528066729  
   NrOfCrOrbitPoints = 10000000  
   period = 1  
   n = 6  
   t = 0.3333332222223279 
   c = -0.1249993954008480 , 0.6495190528380124  
   NrOfCrOrbitPoints = 100000000  
   period = 1  
   n = 7  
   t = 0.3333333222222233 
   c = -0.1249999395400276 , 0.6495190528383258  
   NrOfCrOrbitPoints = 1000000000  
   period = 1  
   n = 8  
   t = 0.3333333322222222 
   c = -0.1249999939540022 , 0.6495190528383290  
   NrOfCrOrbitPoints = 1000000000  
   period = 1  
   n = 9  
   t = 0.3333333332222223 
   c = -0.1249999993954002 , 0.6495190528383290  
   NrOfCrOrbitPoints = 10000000000  
   period = 1  
   n = 10  
   t = 0.3333333333222222 
   c = -0.1249999999395400 , 0.6495190528383290  
   NrOfCrOrbitPoints = 100000000000  
   period = 1  
   n = 11  
   t = 0.3333333333322222 
   c = -0.1249999999939540 , 0.6495190528383290  
   NrOfCrOrbitPoints = 100000  

Please note that all values are converted to double precision. In info procedure :

 printf("   t = %.16f \n", get_double(tt(n)));

Faq

  • Q1. "the Siegel disk looks too smooth near the exterior when the interior digits are deep (starting from n=6 in the cas of 1/3 for instance): maybe the rotation number does not have enough decimals?" Arnaud Cheritat
    • A1 Yes. As n grows time of making image also grows. For last imageit was one day ( next would be probably a week or more) Even then there are not so many points "near the exterior" So I manually join points to get closed curve. Because there is not so many points there then the shape is not precise.  Near fixed points ( inside) there are many points so here image should be precise

references

  1. InfoldingSiegelDisk.gif by Arnaud Chéritat

Captions

Add a one-line explanation of what this file represents

Items portrayed in this file

depicts

16 December 2016

File history

Click on a date/time to view the file as it appeared at that time.

Date/TimeThumbnailDimensionsUserComment
current18:49, 4 January 2017Thumbnail for version as of 18:49, 4 January 2017999 × 999 (55 KB)Soul windsurferbetter quality : more points in outside parts of the image because of inverse iteration. Smaller number of frames . Frames which are hard to draw precisely arre removed
09:42, 17 December 2016Thumbnail for version as of 09:42, 17 December 2016999 × 999 (68 KB)Soul windsurfer0 frame
21:12, 16 December 2016Thumbnail for version as of 21:12, 16 December 2016999 × 999 (64 KB)Soul windsurferUser created page with UploadWizard

The following page uses this file:

Global file usage

The following other wikis use this file:

Metadata