この作品の著作権者である私は、この作品を以下のライセンスで提供します。

このファイルはクリエイティブ・コモンズ 表示-継承 3.0 非移植ライセンスのもとに利用を許諾されています。

あなたは以下の条件に従う場合に限り、自由に

• 共有 – 本作品を複製、頒布、展示、実演できます。
• 再構成 – 二次的著作物を作成できます。

あなたの従うべき条件は以下の通りです。

• 表示 – あなたは適切なクレジットを表示し、ライセンスへのリンクを提供し、変更があったらその旨を示さなければなりません。これらは合理的であればどのような方法で行っても構いませんが、許諾者があなたやあなたの利用行為を支持していると示唆するような方法は除きます。
• 継承 – もしあなたがこの作品をリミックスしたり、改変したり、加工した場合には、あなたはあなたの貢献部分を元の作品とこれと同一または互換性があるライセンスの下に頒布しなければなりません。

https://creativecommons.org/licenses/by-sa/3.0CC BY-SA 3.0 Creative Commons Attribution-Share Alike 3.0 truetrue

• 
• 
  critical orbit, period 7 orbit, center of Siegel disc and ...

rotation number t is :^[1]

${\displaystyle t=[3,2,1000,1,...]=[0;3,2,1000,1\dots ]={\cfrac {1}{3+{\cfrac {1}{2+{\cfrac {1}{1000+{\cfrac {1}{1+\ddots }}}}}}}}}$

In Maxima CAS one can compute it :

kill(all)$
a: [0,3,2,1000,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]$
t:float(cfdisrep(a))$
l:%e^(2*%pi*%i*t)$
c:(l*(1-l/2))/2$
float(rectform(c))

So

t = .2857346725405882
c = 0.113891513213121  +0.595978335936124 i

How one can find limit of [3,2,1000,1,...] ?

Here is explanation of Bill Wood

"I don't know if Maxima knows much about the algebra of continued fractions, but it can be of some help hacking out the manipulation details of a derivation. A most useful fact is that

   [a1, a2, a3, ...] = a1 + 1/[a2, a3, ...]

provided the continued fraction converges. If we apply that three times by hand to [3, 2, 1000, 1, 1, ...] we obtain

   3 + 1/(2 + 1/(1000 + 1/[1, 1, ...]))

Now it is known that [1, 1, ...] converges to the Golden Ratio = (1+sqrt(5))/2. So now we can use Maxima as follows:

 (%i20) 3+(1/(2+(1/(1000+(1/((1+sqrt(5))/2))))));
                                    1
 (%o20)                    ---------------------- + 3
                                  1
                          ------------------ + 2
                               2
                          ----------- + 1000
                          sqrt(5) + 1
 (%i21) factor(%o13);
                              7003 sqrt(5) + 7017
 (%o21)                        -------------------
                              2001 sqrt(5) + 2005
 (%i22) %o21,numer;
 (%o22)                         3.499750279196346

You set a to [0, 3, 2, 1000, 1, 1, ...], which by our useful fact must be the reciprocal of [3, 2, 1000, 1, 1, ...], and indeed the reciprocal of 3.499750279196346 is 0.2857346725405882, which is what your float(t) evaluates to, so we seem to get consistent results.

If all of the continued fractions for the rotation numbers exhibited on the link you provided do end up repeating 1 forever then the method I used above can be used to determine their limits as ratios of linear expressions in sqrt(5)." Bill Wood

"In this example there are 7 big spirals around the Julia set. I guess this is because the rotation number is approximately [3, 2, 0, 0 , 0, ....] = 2/7.

I have checked that there is a 7-periodic point within each spiral, which will be slowly repelling. This shall be the center.

Other spirals will have preperiodic centers with period 7 and any preperiod. Maybe there are other smaller spirals at the Siegel disk belonging to higher periods, which correspond to truncating the continued fraction expansion later. " (Wolf Jung )

In program mandel^[2]

• press key b
• input rotation number .2857346725405882
• switch to dynamic plane ( F2 key)
• increase iteration number ( key n ) up to 65 000
• redraw
• draw critical orbit ( key o) ( to draw full critical orbit, "you must recompile Mandel with a larger value of plotiter in calling drawOrbit() in qmnShell" )
• find periodic point ( key x) for period 7
• make forward iteraion of this point ( key f) to see that it is a periodic point

${\displaystyle z_{0}=-0.291914002678583+0.478643728921302i}$
${\displaystyle z_{1}=-0.029994521062735+0.316532722403284i}$
${\displaystyle z_{2}=0.014598220154870+0.576989841117783i}$
${\displaystyle z_{3}=-0.218812655508313+0.612824385391645i}$
${\displaystyle z_{4}=-0.213783235906926+0.327790873680532i}$
${\displaystyle z_{5}=0.052147928299711+0.455825948583758i}$
${\displaystyle z_{6}=-0.091166375763210+0.643519093703911i}$
${\displaystyle z_{7}=z_{0}=-0.291914002678583+0.478643728921302i}$

Angles of rays landing on period 7 orbit :

9/127
18/127
36/127
72/127
144-127=17/127
34/127
68/127

/* 
   c console  program:
   1. draws Julia setfor Fc(z)=z*z +c
   using DEM/J algorithm ( Distance Esthimation Method for Julia set )

   

   http://www.warwick.ac.uk/~masiay/Research/Siegel.html
   -------------------------------         
   2. technic of creating ppm file is  based on the code of Claudio Rocchini
   http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
   create 24 bit color graphic file ,  portable pixmap file = PPM 
   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
   comments : Adam Majewski 

   gcc s.c -lm -Wall -march=native -O2

   it creates a.out file. Then run it :

   ./a.out

file  s10000.ppm saved

 time ./a.out

real	484m43.557s
user	484m10.020s
sys	0m10.110s

Using Image Magic :
convert s10000.ppm -resize 1000x1000 s10000.png

*/

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

int GiveEscapeTimeJ(double _Zx0, double _Zy0,double C_x, double C_y, int iMax, double _ER2)
{ 
  int i;
  double Zx, Zy;
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
 
  Zx=_Zx0; /* initial value of orbit  */
  Zy=_Zy0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
 
  for (i=0;i<iMax && ((Zx2+Zy2)<_ER2);i++)
    {
      Zy=2*Zx*Zy + C_y;
      Zx=Zx2-Zy2 +C_x;
      Zx2=Zx*Zx;
      Zy2=Zy*Zy;
    };
  return i;
}

 

/*
  estimates distance from point c to nearest point in Julia  set 
  for Fc(z)= z*z + c
  z(n+1) = Fc(zn)  
  this function is based on function  mndlbrot::dist  from  mndlbrot.cpp
  from program mandel by Wolf Jung (GNU GPL )
  http://www.mndynamics.com/indexp.html 

  Hyunsuk Kim  : 
  For Julia sets, z is the variable and c is a constant. Therefore df[n+1](z)/dz = 2*f[n]*f'[n] -- you don't add 1.

  For the Mandelbrot set on the parameter plane, you start at z=0 and c becomes the variable. df[n+1](c)/dc = 2*f[n]*f'[n] + 1. 

*/
double jdist(double Zx, double Zy, double Cx, double Cy ,  int iter_max)
{ 
  int i;
  double x = Zx, /* Z = x+y*i */
    y = Zy, 
    /* Zp = xp+yp*1 = 1  */
    xp = 1, 
    yp = 0, 
    /* temporary */
    nz,  
    nzp,
    /* a = abs(z) */
    a; 
  for (i = 1; i <= iter_max; i++)
    { /* first derivative   zp = 2*z*zp  = xp + yp*i; */
      nz = 2*(x*xp - y*yp) ; 
      yp = 2*(x*yp + y*xp); 
      xp = nz;
      /* z = z*z + c = x+y*i */
      nz = x*x - y*y + Cx; 
      y = 2*x*y + Cy; 
      x = nz; 
      /* */
      nz = x*x + y*y; 
      nzp = xp*xp + yp*yp;
      if (nzp > 1e60 || nz > 1e60) break;
    }
  a=sqrt(nz);
  /* distance = 2 * |Zn| * log|Zn| / |dZn| */
  return 2* a*log(a)/sqrt(nzp); 
}
/* ------------------------------------------------------*/
int main(void)
{ 
  const double Cx=0.113891513213121;
  const double Cy=0.595978335936124;
  /* screen ( integer) coordinate */
  int iX,iY;
  const int iXmax = 1000; 
  const int iYmax = 1000;
  /* world ( double) coordinate = parameter plane*/
  const double ZxMin=-1.5;
  const double ZxMax= 1.5;
  const double ZyMin=-1.5;
  const double ZyMax= 1.5;
  /* */
  double PixelWidth=(ZxMax-ZxMin)/iXmax;
  double PixelHeight=(ZyMax-ZyMin)/iYmax;
  /* color component ( R or G or B) is coded from 0 to 255 */
  /* it is 24 bit color RGB file */
  const int MaxColorComponentValue=255; 
  
  FILE * fp;
  char *filename="s1000.ppm";
  char *comment="# ";/* comment should start with # */

  static unsigned char color[3];

  double Zx0, Zy0;  /* Z0 = Zx0 + Zy0*i */
    
  /*  */
  int LastIteration;
  const int IterationMax=100000;
  /* bail-out value , radius of circle ;  */
  const int EscapeRadius=1/2 + sqrt(1/4+ sqrt(Cx*Cx + Cy*Cy));
  int ER2=EscapeRadius*EscapeRadius;
  double distanceMax=PixelWidth/10; /*jdist( 0,0,Cx,Cy, IterationMax);*/
 

 /*create new file,give it a name and open it in binary mode  */
  fp= fopen(filename,"wb"); /* b -  binary mode */
  /*write ASCII header to the file*/
  fprintf(fp,"P6\n %s\n %d\n %d\n %d\n",comment,iXmax,iYmax,MaxColorComponentValue);

  /* compute and write image data bytes to the file*/
  for(iY=0;iY<iYmax;++iY)
    {
      Zy0=ZyMax - iY*PixelHeight; /* reverse Y  axis */
      if (fabs(Zy0)<PixelHeight/2) Zy0=0.0; /*  */   
      printf("iY/iYmax =   %d / %d \n", iY, iYmax); 
      for(iX=0;iX<iXmax;++iX)
	{    /* initial value of orbit Z0 */
	  Zx0=ZxMin + iX*PixelWidth;
          if (Zx0*Zx0+Zy0*Zy0>ER2) LastIteration=0;
            else LastIteration = GiveEscapeTimeJ(Zx0, Zy0, Cx, Cy, 100, ER2);
            
	  /* compute  pixel color (24 bit = 3 bytes) */
	  if (LastIteration==IterationMax)
	    { /*  interior of Julia set   */
	      color[0]=0;/* this part of code does not work, I can't find the error */
	      color[1]=255;
	      color[2]=0;                           
	    }
	  else /* exterior of Filled-in Julia set  =  */
	    {  double distance=jdist(Zx0,Zy0,Cx,Cy,IterationMax);
	      if (distance<distanceMax)
		{ /*  Julia set  = black */
		  color[0]=0; /* Red*/
		  color[1]=0;  /* Green */ 
		  color[2]=0;/* Blue */
		}
	      else 
		{ /*  exterior of Julia set  = white */
		  color[0]=255;
		  color[1]=255;
		  color[2]=255;                           
		};
	    }
	  /* check the orientation of Z-plane */
	  /* mark first quadrant of cartesian plane*/     
	  /*     if (Z0x>0 && Z0y>0) color[0]=255-color[0];  */
	  /*write color to the file*/
	  fwrite(color,1,3,fp);
	}
    }
  fclose(fp);
  printf("file  %s saved\n", filename);
  
  return 0;
}

1. ↑ Some examples of quadratic polynomial Siegel disks by Davoud Cheraghi
2. ↑ Program mandel by Wolf Jung