File:Parabolic Julia set for internal angle 1 over 15.png

From formulasearchengine
Jump to navigation Jump to search

Original file(1,000 × 1,000 pixels, file size: 105 KB, MIME type: image/png)

This file is from Wikimedia Commons and may be used by other projects. The description on its file description page there is shown below.

Summary

Description
English: Numerical approximation of parabolic Julia set for internal angle 1 over 15. For more detail see long description below.
Date
Source Own program which uses the code by Wolf Jung http://www.mndynamics.com/indexp.html.
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 3.0 Unported 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.

Long description

Discrete dynamical system :

where function f is a complex quadratic polynomial :

Image shows part of dynamical plane with corners :

/* world ( double) coordinate = dynamic plane */
const double ZxMin=-1.5;
const double ZxMax=1.5;
const double ZyMin=-1.5;
const double ZyMax=1.5;

Image shows only Julia set. Note that near alfa fixed point there are 15 components which have a common root point ( z=alfa ). There is also Leau-Fatou flower with 15 petals and center at the alfa fixed point ( it is not draw here but this describes the local dynamic )

Algorithm

This c program can draw parabolic Julia sets for internal angle ( rotation number) from 1/2 to 1/15. One has only to change value of iPeriodChild.

Check every point of dynamic z-plane and compute its color. Points are divided into 2 sets :

  • exterior ( if point escapes , standard bailout test )
  • interior ( if points fall into target set ( trap ) around fixed point alfa.

Setup

Rays landing on alfa fixed points for periods 2 - 40

Compute internal angle ( rotation number) = angle of internal ray of period 1 component of Mandelbrot set :[1]

InternalAngle = 1.0/((double) iPeriodChild); // = 1/15

Compute parameter c :[2]

c = GiveC(InternalAngle, 1.0, 1) ; //

Compute alfa fixed point :[3]

alfa = GiveAlfaFixedPoint(c);

Compute trap for points of interior of Julia set. All points of interior of filled Julia set fall into alfa fixed point thru 15 ( = iPeriodChild) petals. Julia set near alfa fixed point can be aproximated by external rays that land on alfa fixed point. In RaysArray there are description of these rays.

// array of turns
FillRaysArray(iterMax);
FillGapsInRaysArray(); //  many cells of array are empty so fill it using straight lines
OrderRaysArray(); // order rays using their external angles : from smallest to largest. It will make computations easier

Angle of the first external ray that land on alfa fixed point is[4]:

t0 = 1.0/( pow(2.0,iPeriodChild) -1.0); // = 1/32767

One can compute next angles of external rays :

 for (j = 0; j < iPeriodChild ; j++)
  { t *= 2; // t = 2*t
    if (t > 1) t--; // t = t modulo 1 
  }

So numerators are ( 15 numbers) :

1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384

Then compute all rays using backward iteration[5] The good point is that this procedure needs only one angle for computing all 15 rays ( code and algorithm by Wolf Jung)

FillArraySymmetric

For each pixel of dynamic z-plane compute its color using function :

unsigned char GiveColor(unsigned int ix, unsigned int iy)

If points escapes under iteration :

// bailout test
if (Zx2 + Zy2 > ER2) return iExterior; // if escaping stop iteration

give him color of exterior. If not escapes give him color of interior :

Check if it is outside target set around alfa every 15 iteration !!!

Main parameters of the program

  • iPeriodChild ( c is a root point between period 1 and period = iPeriodChild components of Mandelbrot set) It means that c is

from elephant valley

  • iterMax; number here show how smooth will be a boundary ( Julia set ) near parabolic fixed points

Of course more smooth means more time to compute it,

  • iMaxDistance2Alfa = radius of circle around alfa. This circle is a target set ( trap) for points that fall

into alfa fixed point = points of interior of Filled Julia set ( more radius = faster but also maybe distorted image, check it )

Text output of the program

 memory is OK 
i= 0 color = 235  
i= 1 color = 225  
i= 2 color = 215  
i= 3 color = 205  
i= 4 color = 195  
i= 5 color = 185  
i= 6 color = 175  
i= 7 color = 165  
i= 8 color = 155  
i= 9 color = 145  
i= 10 color = 135  
i= 11 color = 125  
i= 12 color = 115  
i= 13 color = 105  
i= 14 color = 95  
landing point of ray for angle = 0.000031 is = (0.360144 ;  0.190271 ) ; iDistnace = 65 
landing point of ray for angle = 0.000061 is = (0.382990 ;  0.154632 ) ; iDistnace = 58 
landing point of ray for angle = 0.000122 is = (0.412261 ;  0.136027 ) ; iDistnace = 53 
landing point of ray for angle = 0.000244 is = (0.440945 ;  0.129740 ) ; iDistnace = 50 
landing point of ray for angle = 0.000488 is = (0.467091 ;  0.131998 ) ; iDistnace = 48 
landing point of ray for angle = 0.000977 is = (0.490240 ;  0.140892 ) ; iDistnace = 47 
landing point of ray for angle = 0.001953 is = (0.509975 ;  0.155724 ) ; iDistnace = 47 
landing point of ray for angle = 0.003906 is = (0.525314 ;  0.176413 ) ; iDistnace = 49 
landing point of ray for angle = 0.007813 is = (0.534324 ;  0.202927 ) ; iDistnace = 51 
landing point of ray for angle = 0.015625 is = (0.533813 ;  0.234439 ) ; iDistnace = 55 
landing point of ray for angle = 0.031251 is = (0.519484 ;  0.267875 ) ; iDistnace = 59 
landing point of ray for angle = 0.062502 is = (0.487597 ;  0.295896 ) ; iDistnace = 65 
landing point of ray for angle = 0.125004 is = (0.439686 ;  0.306138 ) ; iDistnace = 69 
landing point of ray for angle = 0.250008 is = (0.389093 ;  0.286791 ) ; iDistnace = 71 
landing point of ray for angle = 0.500015 is = (0.358634 ;  0.240759 ) ; iDistnace = 70 
landing point of ray for angle = 1.000031 is = (0.360144 ;  0.190271 ) ; iDistnace = 65 
in ray j= 0 there is a gap from i= 0 to i= 65 iStep = 65 ; 
in ray j= 0 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 1 there is a gap from i= 0 to i= 58 iStep = 58 ; 
in ray j= 1 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 2 there is a gap from i= 0 to i= 53 iStep = 53 ; 
in ray j= 2 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 3 there is a gap from i= 0 to i= 50 iStep = 50 ; 
in ray j= 3 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 4 there is a gap from i= 0 to i= 48 iStep = 48 ; 
in ray j= 4 there is a gap from i= 110 to i= 112 iStep = 2 ; dStep = 0.000471
in ray j= 4 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 5 there is a gap from i= 0 to i= 47 iStep = 47 ; 
in ray j= 5 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 6 there is a gap from i= 0 to i= 47 iStep = 47 ; 
in ray j= 6 there is a gap from i= 113 to i= 115 iStep = 2 ; dStep = 0.000368
in ray j= 6 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 7 there is a gap from i= 0 to i= 49 iStep = 49 ; 
in ray j= 7 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 8 there is a gap from i= 0 to i= 51 iStep = 51 ; 
in ray j= 8 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 9 there is a gap from i= 0 to i= 55 iStep = 55 ; 
in ray j= 9 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 10 there is a gap from i= 0 to i= 59 iStep = 59 ; 
in ray j= 10 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 11 there is a gap from i= 0 to i= 65 iStep = 65 ; 
in ray j= 11 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 12 there is a gap from i= 0 to i= 69 iStep = 69 ; 
in ray j= 12 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 13 there is a gap from i= 0 to i= 71 iStep = 71 ; 
in ray j= 13 there is a gap from i= 118 to i= 120 iStep = 2 ; 
in ray j= 14 there is a gap from i= 0 to i= 70 iStep = 70 ; 
in ray j= 14 there is a gap from i= 118 to i= 120 iStep = 2 ; 
axis of symmetry 
File 100000000.000000.pgm saved. 
 find boundaries in data array using  Sobel filter
File 100000000.001000.pgm saved. 
copy boundaries from edge array to data array 
File 100000000.002000.pgm saved. 
File 100000000.003000.pgm saved. 
InternalAngle  = 0.066667 
Cx  = 0.289490 
Cy  = 0.017582 
alfax  = 0.456773 
alfay  = 0.203368 
iHeight  = 2001 
PixelWidth  = 0.001500 
distorsion of image  = 1.000000 
iterMax = 100000000 
dMaxDistance2Alfa= 0.180000
iMaxDistance2Alfa= 120

real	23m42.638s
user	31m19.890s
sys	0m4.040s

C code

/*

  Adam Majewski
  fraktal.republika.pl

  c console progam using 
  * symmetry
  * openMP

  draw  parabolic Julia set 

  and saves it to pgm files ( different versions )

 

  gcc r7.c -lm -Wall -fopenmp -march=native 
  time ./a.out

  main parameters :
  - iPeriodChild ( c is a root point between period 1 and period = iPeriodChild components of Mandelbrot set) It means that c is
  from elephant valley

  - FillRaysArray(10000000); number here show how smooth will be a boundary ( Julia set ) near parabolic fixed points
  Of course more smooth means more time to compute iter

  -  iMaxDistance2Alfa = radius of circle around alfa. This circle is a target set ( trap) for points that fall 
  into alfa fixed point  = points of interior of Filled Julia set ( more radius = faster but also maybe distorted image, check it )

*/

#include <stdio.h>
#include <stdlib.h> // malloc
#include <string.h> // strcat
#include <math.h> // M_PI; needs -lm also 
#include <complex.h>
#include <omp.h> // OpenMP; needs also -fopenmp

/* --------------------------------- global variables and constans ------------------------------------------------------------ */
// iPeriodChild of secondary component joined by root point
#define iPeriodChild 15

//  unsigned int denominator;  denominator = iPeriodChild;
double InternalAngle;
unsigned char Colors[iPeriodChild]; //={255,230,180, 160,140,120,100}; // NumberOfPetal of colors = iPeriodChild
unsigned char iExterior = 245;
unsigned char iPetal = 255;

// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1 
unsigned int ix, iy; // var
unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
unsigned int ixMax ; //
unsigned int iWidth ; // horizontal dimension of array
unsigned int ixAxisOfSymmetry  ; // 
unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
unsigned int iyMax ; //
unsigned int iyAxisOfSymmetry  ; // 
unsigned int iyAbove ; // var, measured from 1 to (iyAboveAxisLength -1)
unsigned int iyAboveMin = 1 ; //
unsigned int iyAboveMax ; //
unsigned int iyAboveAxisLength ; //
unsigned int iyBelowAxisLength ; //
unsigned int iHeight = 2000; //  odd NumberOfPetal !!!!!! = (iyMax -iyMin + 1) = iyAboveAxisLength + iyBelowAxisLength +1
// The size of array has to be a positive constant integer 
unsigned int iSize ; // = iWidth*iHeight; 

// memmory 1D arrays 
unsigned char *data;
unsigned char *edge;

// unsigned int i; // var = index of 1D array
unsigned int iMin = 0; // Indexes of array starts from 0 not 1
unsigned int iMax ; // = i2Dsize-1  = 
// The size of array has to be a positive constant integer 
// unsigned int i1Dsize ; // = i2Dsize  = (iMax -iMin + 1) =  ;  1D array with the same size as 2D array

/* world ( double) coordinate = dynamic 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;
 
double ratio ;

// complex numbers of parametr plane 
double Cx; // c =Cx +Cy * i
double Cy;
double complex c; // 

double complex alfa; // alfa fixed point alfa=f(alfa)
double alfax,alfay;

unsigned long int iterMax  =100000000; //iHeight*100;
// target set for escaping points is a exterior of circle with center in origin
double ER = 2.0; // Escape Radius for bailout test 
double ER2;

#define iMaxDistance2Alfa  120 // distance point to alfa fixed point in pixels where PixelWidth = 0.003; // 3 world units / 1000 pixels 
double dMaxDistance2Alfa2; // = (iMaxDistance2Alfa*PixelWidth)^2

//
// target set for points falling into alfa fixed point
// is a circle around alfa fixed point
// with radius = AR
//double AR ; // radius of target set around alfa fixed point in world coordinate AR = PixelWidth*TargetWidth; 
//double AR2; // =AR*AR;

double TwoPi=2*M_PI;

// array with angles in turns of points of periodic rays landing on alfa fixed point :
// contains  iCrDistance x iPeriodChild  angles 
double RaysTurns[iMaxDistance2Alfa][iPeriodChild];

/* ------------------------------------------ functions -------------------------------------------------------------*/

// gives argument of complex number in turns 
// realted with alfa fixed point
// http://en.wikipedia.org/wiki/Turn_%28geometry%29
double GiveTurn(double zx, double zy)
{
  double argument;

  argument =atan2(zy-alfay, zx-alfax);// carg(zx-alfax +(zy-alfay)*I ); // alfa !!!;   argument in radians from -pi to pi
  if (argument<0) argument=argument + TwoPi; //   argument in radians from 0 to 2*pi
  return argument/TwoPi ; // argument in turns from 0.0 to 1.0
}

/* 
   principal square  root of complex number 
   http://en.wikipedia.org/wiki/Square_root

   z1= I;
   z2 = root(z1);
   printf("zx  = %f \n", creal(z2));
   printf("zy  = %f \n", cimag(z2));
*/
double complex root(double complex z)
{ 
  double x = creal(z);
  double y = cimag(z);
  double u;
  double v;
  double r = sqrt(x*x + y*y); 
  
  v = sqrt(0.5*(r - x));
  if (y < 0) v = -v; 
  u = sqrt(0.5*(r + x));
  return u + v*I;
}

double complex preimage(double complex z1, double complex z2,  double complex c)
{ 
  double complex zPrev;

  zPrev = root(creal(z1) - creal(c) + ( cimag(z1) - cimag(c))*I);
  // choose one of 2 roots 
  if (creal(zPrev)*creal(z2) + cimag(zPrev)*cimag(z2) > 0) 
    return zPrev ;     //  u+v*i
  else return -zPrev; // -u-v*i
}

// This function  works for periodic angles.
// You must determine the period n before calling this function.
// draws all "period" external rays 

// based on the code for backward iteration for drawing external ray see QmnPlane::backRay()
// by Wolf Jung  http://www.mndynamics.com/indexp.html

double complex FillRaysArray(int IterMax  )
{
  
  double xend ; // re of the endpoint of the ray
  double yend; // im of the endpoint of the ray
  const double R = 10000; // very big radius = near infinity
  int j; // number of ray 
  int iter; // index of backward iteration
  double t,t0;
  
  
  double complex zPrev;
  double u,v; // zPrev = u+v*I
  double complex zNext;

  int iDistance ; // dDistance/PixelWidth
  
  // fill array with negative number 
  for (iDistance=0; iDistance<iMaxDistance2Alfa; ++iDistance)
    for (j=0; j<iPeriodChild; ++j)
      RaysTurns[iDistance][j]=-1.0;

  t0 = 1.0/( pow(2.0,iPeriodChild) -1.0); // http://fraktal.republika.pl/mset_external_ray_m.html
  t=t0;

  /* dynamic 1D arrays for coordinates ( x, y) of points with the same R on preperiodic and periodic rays  */
  double *RayXs, *RayYs;
  int iLength = iPeriodChild+2; // length of arrays ?? why +2
  //  creates arrays :  RayXs and RayYs  and checks if it was done
  RayXs = malloc( iLength * sizeof(double) );
  RayYs = malloc( iLength * sizeof(double) );
  if (RayXs == NULL || RayYs==NULL)
    {
      fprintf(stderr,"Could not allocate memory");
      getchar(); 
      return 1; // error
    }
  

  //  starting points on preperiodic and periodic rays 
  //  with angles t, 2t, 4t...  and the same radius R
  for (j = 0; j < iPeriodChild ; j++)
    { // z= R*exp(2*Pi*t)
      RayXs[j] = R*cos((2*M_PI)*t); 
      RayYs[j] = R*sin((2*M_PI)*t);

    
      t *= 2; // t = 2*t
      if (t > 1) t--; // t = t modulo 1 
    }
  zNext = RayXs[0] + RayYs[0] *I;

  // ???
  // z[k] is n-periodic. So it can be defined here explicitly as well.
  RayXs[iPeriodChild] = RayXs[0]; 
  RayYs[iPeriodChild] = RayYs[0];
  

  //   backward iteration of each point z
  for (iter = -10; iter <= IterMax; iter++)
    { 
     	
      for (j = 0; j < iPeriodChild; j++) // period +preperiod
	{ // u+v*i = sqrt(z-c)   backward iteration in fc plane 
	  zPrev = root(RayXs[j+1] - creal(c)+(RayYs[j+1] - cimag(c))*I ); // , u, v
	  u=creal(zPrev);
	  v=cimag(zPrev);
                
	  // choose one of 2 roots: u+v*i or -u-v*i
	  if (u*RayXs[j] + v*RayYs[j] > 0) 
	    { RayXs[j] = u; RayYs[j] = v; } // u+v*i
	  else { RayXs[j] = -u; RayYs[j] = -v; } // -u-v*i

	  //if inside trap !! save turns to the array
	  iDistance = (int)(sqrt((RayXs[j]-alfax)*(RayXs[j]-alfax) +  (RayYs[j]-alfay)*(RayYs[j]-alfay))/PixelWidth);
	  if ( iDistance < iMaxDistance2Alfa )
                   
	    {
                   
	      RaysTurns[iDistance][j]= GiveTurn(  RayXs[j],  RayYs[j]);
	    }
                  
                 
                
	} // for j ...

      // ???
      // z[k] is n-periodic. So it can be defined here explicitly as well.
      RayXs[iPeriodChild] = RayXs[0]; 
      RayYs[iPeriodChild] = RayYs[0];   
          
     
    }

  // check 
  t=t0;
  for (j = 0; j < iPeriodChild + 1; j++)
    {
      // aproximate end of ray by straight line to its landing point here = alfa fixed point
      //dDrawLine(RayXs[j],RayYs[j], creal(alfa), cimag(alfa), 0, data); 
      iDistance = (int)(sqrt((RayXs[j]-alfax)*(RayXs[j]-alfax) +  (RayYs[j]-alfay)*(RayYs[j]-alfay))/PixelWidth);
      printf("landing point of ray for angle = %f is = (%f ;  %f ) ; iDistnace = %d \n",t, RayXs[j], RayYs[j],  iDistance);
      t *= 2; // t = 2*t
    } // end of the check 

  // last point of a ray 0
  xend = RayXs[0];
  yend = RayYs[0];

  

  // free memmory
  free(RayXs);
  free(RayYs);

  return  xend + yend*I; // return last point or ray for angle t 
}

// aproximate not computed ( negative values)
int FillGapsInRaysArray()

{
  // indexes of the array 
  int i; // index of the pixel on the ray
  int iMini, iMaxi; // gap border
  int j; // index of the ray
 
  double dStep;
  int iStep;

  // fill empty values = change negative with positive !!!!!!!!!!
  // negative value = not filled = gaps
  // positive value = filled ( computed )

  // check  every ray 
  for (j = 0; j < iPeriodChild ; j++)
    {

      // start from 0 wher are negative values 
      // because of slow dynamics 
      i=0;
      // find first positive value 
      while (RaysTurns[i][j]<0.0) i+=1;
      iMaxi = i;// here positive

      // copy first positive value  to all cells before it = straight line 
      for (i=0; i<iMaxi; i++) RaysTurns[i][j] = RaysTurns[iMaxi][j];
      // go back 
      i = iMaxi; // positive
      printf("in ray j= %d there is a gap from i= 0 to i= %d iStep = %d ; \n",j,  iMaxi, iMaxi);
      // rest of the array
      do {
	// go thru all positive 
	while (RaysTurns[i][j]>0.0 && i<iMaxDistance2Alfa-1) i+=1;
	iMini=i-1; // here is positive, lower border of the gap 
	// here value is negative : RaysTurns[i][j]<0.0
	do i+=1; while (RaysTurns[i][j]<0.0 && i<iMaxDistance2Alfa-1  ); 
	iMaxi= i; // positive , upper border of the gap
	iStep= iMaxi-iMini;
	// step between RaysTurns[iMini][j] and RaysTurns[iMaxi][j]
	if (RaysTurns[iMaxi][j]>RaysTurns[iMini][j] ) 
	  dStep = (RaysTurns[iMaxi][j]-RaysTurns[iMini][j])/iStep;
	else dStep = (RaysTurns[iMini][j]-RaysTurns[iMaxi][j])/iStep;
	// step between values 
	if (iMaxi<iMaxDistance2Alfa-1 && RaysTurns[iMaxi-1][j]< 0.0) // array is numberd from 0 to ...
	  {
	    printf("in ray j= %d there is a gap from i= %d to i= %d iStep = %d ; dStep = %f\n",j, iMini, iMaxi, iStep, dStep );
	    // fill the gap 
	    i= iMini;
	    // while (i<iMaxi-1) { i+=1; RaysTurns[i][j]=RaysTurns[iMini][j] + (i-iMini)*dStep;} 
	    i=iMaxi;
	  }
        else {
	  printf("in ray j= %d there is a gap from i= %d to i= %d iStep = %d ; \n",j,  iMini, iMaxi, iStep);
	  // copy last positive value  to all cells after it = straight line 
	  for (i=iMini+1; i<iMaxi; i++) RaysTurns[i][j] = RaysTurns[iMini][j];
	} 
       
      } while (i<iMaxDistance2Alfa);
  
    } // j

  return 0;
}

// row of array RaysTurns [j] 
// will be ordered from :
// lowest angle  RaysTurns[j][0]>0.0
// to maximal angle RaysTurns[j][iMaxDistance2Alfa-1] < 1.0
int OrderRaysArray()

{ 

  double row[iPeriodChild];
  double SmallestAngle;
  int i,j;
  
  int iSmallest;

  // find smallest angle in Turns array 
  SmallestAngle = RaysTurns[0][0];
  for (i=0; i<iPeriodChild; ++i) if (RaysTurns[0][i]<SmallestAngle) {SmallestAngle=RaysTurns[0][i]; iSmallest=i;}

  for (i=0; i<iMaxDistance2Alfa; ++i)  

    {
      for (j = 0; j < iPeriodChild ; j++)
	row[j]= RaysTurns[i][(iSmallest+j) % iPeriodChild]; // copy to row arrays
      for (j = 0; j < iPeriodChild ; j++)
	// copy from row to RayTurns   
	RaysTurns[i][j] = row[j];
   

    }
 

  return 0;
}

int SaveArray2File(double a[iMaxDistance2Alfa][iPeriodChild])
{

  int d,p;

  FILE * fp;
  fp= fopen("file.txt","wb"); /*create new file,give it a name and open it in binary mode  */
  for (d=0; d<iMaxDistance2Alfa ; ++d)
    { fprintf(fp, " i= %d ", d);
      {  for (p=0; p<iPeriodChild ; ++p) fprintf(fp,"turn = %f ", a[d][p]);
	fprintf(fp, "\n");}  } 
  fclose(fp); 

  return 0;

}

/* 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

*/
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int iPeriod)
{
  //0 <= InternalRay<= 1
  //0 <= InternalAngleInTurns <=1
  double t = InternalAngleInTurns *2*M_PI; // from turns to radians
  double R2 = InternalRadius * InternalRadius;
  //double Cx, Cy; /* C = Cx+Cy*i */
  switch ( iPeriod ) // 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
      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 // general form  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)) 
double complex GiveAlfaFixedPoint(double complex c)
{
  double dx, dy; //The discriminant is the  d=b^2- 4ac = dx+dy*i
  double r; // r(d)=sqrt(dx^2 + dy^2)
  double sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i
  double 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;
}

// colors of components interior = shades of gray
int InitColors()
{
  int i;
  int iMax = iPeriodChild; // uses global var iPeriodChild and Colors
  unsigned int iStep;

  iStep=150/iPeriodChild;

  for (i = 1; i <= iMax; ++i)
    {Colors[i-1] = iExterior -i*iStep; 
      printf("i= %d color = %i  \n",i-1, Colors[i-1]);}
  return 0;
}

/* -------------------------------------------------- SETUP --------------------------------------- */
int setup()

{

  
  /* 2D array ranges */
  if (!(iHeight % 2)) iHeight+=1; // it sholud be even NumberOfPetal (variable % 2) or (variable & 1)
  iWidth = iHeight;
  iSize = iWidth*iHeight; // size = NumberOfPetal of points in array 
  
  // iy
  iyMax = iHeight - 1 ; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].
  iyAboveAxisLength = (iHeight -1)/2;
  iyAboveMax = iyAboveAxisLength ; 
  iyBelowAxisLength = iyAboveAxisLength; // the same 
  iyAxisOfSymmetry = iyMin + iyBelowAxisLength ; 
  
  // ix
  ixMax = iWidth - 1;

  /* 1D array ranges */
  // i1Dsize = i2Dsize; // 1D array with the same size as 2D array
  iMax = iSize-1; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].
  
  /* Pixel sizes */
  PixelWidth = (ZxMax-ZxMin)/ixMax; //  ixMax = (iWidth-1)  step between pixels in world coordinate 
  PixelHeight = (ZyMax-ZyMin)/iyMax;
  ratio = ((ZxMax-ZxMin)/(ZyMax-ZyMin))/((float)iWidth/(float)iHeight); // it should be 1.000 ...
   
  // for numerical optimisation 
  ER2 = ER * ER;
  dMaxDistance2Alfa2 = (iMaxDistance2Alfa*PixelWidth)*(iMaxDistance2Alfa*PixelWidth);// AR2 
   
  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc( iSize * sizeof(unsigned char) );
  edge = malloc( iSize * sizeof(unsigned char) );
  if (edge == NULL || edge == NULL )
    {
      fprintf(stderr," Could not allocate memory\n");
      return 1;
    }
  else fprintf(stderr," memory is OK \n");

  
  InitColors();
  

  // 
  InternalAngle = 1.0/((double) iPeriodChild); // 
  c = GiveC(InternalAngle, 1.0, 1) ;
  Cx=creal(c);
  Cy=cimag(c);
  //
  alfa = GiveAlfaFixedPoint(c);
  alfax=creal(alfa);
  alfay=cimag(alfa);

  // array of turns 
  
  FillRaysArray(iterMax);
  //printf(" dist2 alfa = %d \n", (int)(((creal(z)-alfax)*(creal(z)-alfax) +  (cimag(z)-alfay)*(cimag(z)-alfay))/PixelWidth));
  FillGapsInRaysArray(); 
  OrderRaysArray();
  SaveArray2File(RaysTurns);
  
 
  return 0;

}

// from screen to world coordinate ; linear mapping
// uses global cons
double GiveZx(unsigned int ix)
{ return (ZxMin + ix*PixelWidth );}

// uses global cons
double GiveZy(unsigned int iy)
{ return (ZyMax - iy*PixelHeight);} // reverse y axis

// all points of interior fall into parabolic fixed point z=alfa
// thru iPeriodChild petals
// boundaries of petals are aproximated by periodic rays
// landing on the alfa fixed point 
unsigned char GiveColorOfInterior(double x, double y, double distance2alfa2)
 
{
  double angle;
  int iDistance2Alfa;

 
  iDistance2Alfa = (int)(sqrt(distance2alfa2)/PixelWidth); 
 
  int i;
  

  
  //  
  angle=GiveTurn(x,y); // z-alfa !!!!
  if (angle<RaysTurns[iDistance2Alfa][0] || angle>RaysTurns[iDistance2Alfa][iPeriodChild-1])
    return Colors[0]; 
  for(i=1;i<iPeriodChild-1;++i)  if (angle<RaysTurns[iDistance2Alfa][i]) return Colors[i];

  /*j=0;*/
  /*  //printf("i= %d\n",i);*/
  /*   while (angle > RaysTurns[iDistance2Alfa][j] ) */
  /*    j+=1; */
  return Colors[iPeriodChild-1];  
 
    
 
}

unsigned char GiveColor(unsigned int ix, unsigned int iy)
{ 
  // check behavour of z under fc(z)=z^2+c
  // using 2 target set:
  // 1. exterior or circle (center at origin and radius ER ) 
  // as a target set containing infinity = for escaping points ( bailout test)
  // for points of exterior of julia set
  // 2. interior of circle with center = alfa and radius dMaxDistance2Alfa
  // as a target set for points of interior of Julia set 
  double Zx, Zy; //  Z= Zx+ZY*i;
  double Zx2, Zy2;
  int i;
  //int j; // iteration = fc(z)
  double d2 ; /* d2= (distance from z to Alpha)^2   */
  double dX,dY ; // d = dx+dy*I
  
  // from screen to world coordinate 
  Zx = GiveZx(ix);
  Zy = GiveZy(iy);
  /* distance from z to Alpha  */
  dX=Zx-alfax;
  dY=Zy-alfay;
  d2=dX*dX+dY*dY;
  // if inside target set around attractor ( alfa fixed point )
  while (d2>dMaxDistance2Alfa2)
    {
      for(i=0;i<iPeriodChild ;++i) // iMax = period !!!!
	{ 
	  Zx2 = Zx*Zx; 
	  Zy2 = Zy*Zy;
       
	  // bailout test 
	  if (Zx2 + Zy2 > ER2) return iExterior; // if escaping stop iteration
       
	  // if not escaping or not attracting then iterate = check behaviour
	  // new z : Z(n+1) = Zn * Zn  + C
	  Zy = 2*Zx*Zy + Cy; 
	  Zx = Zx2 - Zy2 + Cx; 
	} 
      /* distance from z to Alpha  */
      dX=Zx-alfax;
      dY=Zy-alfay;
      d2=dX*dX+dY*dY;
      // if inside target set around attractor ( alfa fixed point )
      
      
    }
  return GiveColorOfInterior( Zx, Zy, d2);//iPetal; // not escaping and not in attracting target set  , probably never here (:-))
}

/* -----------  array functions -------------- */

/* gives position of 2D point (iX,iY) in 1D array  ; uses also global variable iWidth */
unsigned int Give_i(unsigned int ix, unsigned int iy)
{ return ix + iy*iWidth; }
//  ix = i % iWidth;
//  iy = (i- ix) / iWidth;
//  i  = Give_i(ix, iy);

// plots raster point (ix,iy) 
int PlotPoint(unsigned int ix, unsigned int iy, unsigned char iColor)
{
  unsigned i; /* index of 1D array */
  i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
  data[i] = iColor;

  return 0;
}

// fill array 
// uses global var :  ...
// scanning complex plane 
int FillArray(unsigned char data[] )
{
  unsigned int ix, iy; // pixel coordinate 

  // for all pixels of image 
  for(iy = iyMin; iy<=iyMax; ++iy) 
    { printf(" %d z %d\n", iy, iyMax); //info 
      for(ix= ixMin; ix<=ixMax; ++ix) PlotPoint(ix, iy, GiveColor(ix, iy) ); //  
    } 
   
  return 0;
}

// fill array using symmetry of image 
// uses global var :  ...
int FillArraySymmetric(unsigned char data[] )
{
   
  unsigned char Color; // gray from 0 to 255 

  printf("axis of symmetry \n"); 
  iy = iyAxisOfSymmetry; 
#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)
  for(ix=ixMin;ix<=ixMax;++ix) PlotPoint(ix, iy, GiveColor(ix, iy));

  /*
    The use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.
    The use of ‘private(variable, variable2)’ specifies that these variables should have a seperate instance in each thread.
  */

#pragma omp parallel for schedule(dynamic) private(iyAbove,ix,iy,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)

  // above and below axis 
  for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
    {// printf(" %d from %d\r", iyAbove, iyAboveMax); //info 
      for(ix=ixMin; ix<=ixMax; ++ix) 

	{ // above axis compute color and save it to the array
	  iy = iyAxisOfSymmetry + iyAbove;
	  Color = GiveColor(ix, iy);
	  PlotPoint(ix, iy, Color ); 
	  // below the axis only copy Color the same as above without computing it 
	  PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color ); 
	} 
    }  
  return 0;
}

int AddBoundaries(unsigned char data[])
{

  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
  /* sobel filter */
  unsigned char G, Gh, Gv; 
 
  

  printf(" find boundaries in data array using  Sobel filter\n");   
#pragma omp parallel for schedule(dynamic) private(i,iY,iX,Gv,Gh,G) shared(iyMax,ixMax, ER2)
  for(iY=1;iY<iyMax-1;++iY){ 
    for(iX=1;iX<ixMax-1;++iX){ 
      Gv= data[Give_i(iX-1,iY+1)] + 2*data[Give_i(iX,iY+1)] + data[Give_i(iX-1,iY+1)] - data[Give_i(iX-1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[Give_i(iX+1,iY-1)];
      Gh= data[Give_i(iX+1,iY+1)] + 2*data[Give_i(iX+1,iY)] + data[Give_i(iX-1,iY-1)] - data[Give_i(iX+1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[Give_i(iX-1,iY-1)];
      G = sqrt(Gh*Gh + Gv*Gv);
      i= Give_i(iX,iY); /* compute index of 1D array from indices of 2D array */
      if (G==0) {edge[i]=255;} /* background */
      else {edge[i]=0;}  /* boundary */
    }
  }

  // copy boundaries from edge array to data array 
  //for(iY=1;iY<iyMax-1;++iY){ 
  //   for(iX=1;iX<ixMax-1;++iX){i= Give_i(iX,iY); if (edge[i]==0) data[i]=0;}}

  return 0;
}

// Check Orientation of image : first quadrant in upper right position
// uses global var :  ...
int CheckOrientation(unsigned char data[] )
{
  unsigned int ix, iy; // pixel coordinate 
  double Zx, Zy; //  Z= Zx+ZY*i;
  unsigned i; /* index of 1D array */
  for(iy=iyMin;iy<=iyMax;++iy) 
    {
      Zy = GiveZy(iy);
      for(ix=ixMin;ix<=ixMax;++ix) 
	{

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

	}
    }
   
  return 0;
}

int CopyBoundaries()
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
 
 
  printf("copy boundaries from edge array to data array \n");
  for(iY=1;iY<iyMax-1;++iY)
    for(iX=1;iX<ixMax-1;++iX)
      {i= Give_i(iX,iY); if (edge[i]==0) data[i]=0;}
 
 
 
  return 0;
}

int DrawCriticalOrbit(unsigned int IterMax)
{
 
  unsigned int ix, iy; // pixel coordinate 
  double Zx=0.0;
  double Zy=0.0; //  Z= Zx+ZY*i;
  double Zx2=0.0;
  double Zy2=0.0;
  unsigned int i; /* index of 1D array */
  unsigned int j;

  // draw critical point  
  ix = (int)((Zx-ZxMin)/PixelWidth);
  iy = (int)((ZyMax-Zy)/PixelHeight); // reverse y axis
  i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
  data[i]=255-data[i]; 

  // iterate
  for (j = 1; j <= IterMax; j++) //larg number of iteration s
    {  Zx2 = Zx*Zx; 
      Zy2 = Zy*Zy;
       
      // bailout test 
      if (Zx2 + Zy2 > ER2) return iExterior; // if escaping stop iteration
       
      // if not escaping iterate
      // Z(n+1) = Zn * Zn  + C
      Zy = 2*Zx*Zy + Cy; 
      Zx = Zx2 - Zy2 + Cx;
      //compute integer coordinate  
      ix = (int)((Zx-ZxMin)/PixelWidth);
      iy = (int)((ZyMax-Zy)/PixelHeight); // reverse y axis
      i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
      data[i]=255-data[i];   // mark the trap

    }
  return 0;
}

// save data array to pgm file 
int SaveArray2PGMFile( unsigned char data[], double t)
{
  
  FILE * fp;
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name [100]; /* name of file */
  sprintf(name,"%f", t); /*  */
  char *filename =strcat(name,".pgm");
  char *comment="# ";/* 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,iWidth,iHeight,MaxColorComponentValue);  /*write header to the file*/
  fwrite(data,iSize,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);

  return 0;
}

int info()
{
  // diplay info messages
  printf("InternalAngle  = %f \n", InternalAngle);
  printf("Cx  = %f \n", Cx); 
  printf("Cy  = %f \n", Cy);
  // 
  printf("alfax  = %f \n", creal(alfa));
  printf("alfay  = %f \n", cimag(alfa));
  printf("iHeight  = %d \n", iHeight);
  printf("PixelWidth  = %f \n", PixelWidth);
  
  printf("distorsion of image  = %f \n", ratio);
  printf("iterMax = %lu \n", iterMax);
  printf("dMaxDistance2Alfa= %f\n", sqrt(dMaxDistance2Alfa2));
  printf("iMaxDistance2Alfa= %d\n", iMaxDistance2Alfa);
 

  // ------------------------------------
  return 0;
}

/* -----------------------------------------  main   -------------------------------------------------------------*/
int main()
{

  
  setup();
  
  

  // here are procedures for creating image file
  
  //FillArray( data ); // no symmetry
  FillArraySymmetric(data);
  SaveArray2PGMFile(data , iterMax+0.000); // save edge array (only boundaries) to pgm file 
  AddBoundaries(data);
  //CheckOrientation( data );
  SaveArray2PGMFile(edge ,iterMax+0.001); // save edge array (only boundaries) to pgm file 
  CopyBoundaries();
  SaveArray2PGMFile(data , iterMax+0.002); // save edge array (only boundaries) to pgm file 
 
  DrawCriticalOrbit(1000000);
  SaveArray2PGMFile(data , iterMax+0.003); // save edge array (only boundaries) to pgm file
  //
  free(data);
  free(edge);
  
  //
  info();

  return 0;
}

Image Magic code

convert 100000000.001000.pgm -resize 1000x1000 a.png

References

  1. wikibooks : Fractals/Iterations_in_the_complex_plane/Mandelbrot_set#Internal_angle
  2. wikibooks : Fractals /Iterations_in_the_complex_plane/parabolic#Period_1
  3. wikibooks:Fractals/Mathematics/Period#Fixed_points_.28_period_.3D_1_.29
  4. external ray of root points
  5. wikibooks : Fractals/Iterations_in_the_complex_plane/Julia_set#backwards_iteration

Captions

Add a one-line explanation of what this file represents

Items portrayed in this file

depicts

23 June 2013

File history

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

Date/TimeThumbnailDimensionsUserComment
current20:09, 23 June 2013Thumbnail for version as of 20:09, 23 June 20131,000 × 1,000 (105 KB)wikimediacommons>Soul windsurferUser created page with UploadWizard

There are no pages that use this file.