File:Parabolic Julia set for internal angle 1 over 5 - new method.png

Size of this preview: 600 × 600 pixels. Other resolutions: 240 × 240 pixels | 480 × 480 pixels | 768 × 768 pixels | 1,024 × 1,024 pixels | 1,993 × 1,993 pixels.

Original file (1,993 × 1,993 pixels, file size: 33 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: Parabolic Julia set for fc(z) = z^2 + c where c is on boundary of main cardioid with internal angle 1 over 5. Drawing : new method

Date

26 December 2012

Source

Own work

Author

Adam majewski

Permission
(Reusing this file)

I, the copyright holder of this work, hereby publish it under the following licenses:

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.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of the license is included in the section entitled GNU Free Documentation License.

You may select the license of your choice.

PNG development

This PNG graphic was created with C.

Algorithm

First step is to divide points z of complex dynamical plane into 2 subsets ( using bailout test; see function GiveColor ) :

escaping points = exterior
non-escaping points interior and boundary = filled Julia set

iter  =GiveLastIteration(Zx, Zy );

if (iter< iterMax ) // point escapes   
   { color = 245; } // exterior 
   else // Filled Julia Set

Second step is to divide points of interior into 5 subsets using :

target set
fall-in test

Target set features :

center = alfa.
radius = AR
is divided into 5 subsets ( petals) by 5-arm star
all points of interior fall into fixed point alfa

Fall-in test : iterate z and check when point z is inside target set.

Note that test is done after every $f^{period}$ not after every $f$

 dX=Zx-_ZAx;
 dY=Zy-_ZAy;
 d=dX*dX+dY*dY; // d = square of distance from zn to alfa

 while ( d>_AR2) // check if z is inside target 
    {
      for(i=0;i<period ;++i) // iMax = period !!!!

Find in which subset of target set point z falls using relative angle ( see function GiveNumber ) :

  if ((0.0<angle) && (angle<0.2)) return 0;
  if ((0.2<angle) && (angle<0.4)) return 1;
  if ((0.4<angle) && (angle<0.6)) return 2;
  if ((0.6<angle) && (angle<0.8)) return 3;
  return 4; // angle > 0.8
</source >

The star is turned <ref>Complex Dynamics by  Lennart Carleson,Theodore Williams Gamelin  page 40</ref> so
<syntaxhighlight copy line lang=c>
 angle =GiveTurn(dX +dY*I) - 0.17298227404701; // this turn offset is computed with Maxima CAS file

So range of angle in which point zn falls into target set gives a color of subset of interior of Julia set :

 unsigned char colors[]={255,230,180, 160,140}; // number of colors >= period
 // ....
 int number =GiveNumber(Zx,Zy, Cx, Cy, period, AR2, creal(alfa), cimag(alfa));
 color = colors[number];
</source >
'''Third step''' : when all points of plane are colored apply edge detection with Sobel filter ( see function  AddBoundaries )

'''Forth step''' : save memory array to pgm file

'''Five step''' : convert pgm to png using Image Magic :

 convert a.pgm a.png

== {{int:license-header}} ==
{{self|cc-by-sa-3.0|GFDL}}
==Maxima CAS src code ==
This program computes turn offset = angle in turns. This angle equal to rotation of 5-arm star with center in alfa fixed point. 

It works here but for other values ( period , q) do not. Help is welcome to improve this code !!!!

<syntaxhighlight copy line lang=wikitext>
kill(all);
remvalue(all);



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

f(z,c):=z*z+c;

fn(p, z, c) :=
  if p=0 then z
  elseif p=1 then f(z,c)
  else f(fn(p-1, z, c),c);


/* gives parameter c for complex quadratic polynomial */
GiveC(angle,radius):=
(
 [w,c],
 w:radius*%e^(%i*angle*2*%pi),  /* point of circle */
 /* conformal map  from circle to cardioid ( boundary of period 1 component of Mandelbrot set */
 c:w/2-w*w/4, 
 float(rectform(c))    /* point on boundary of period 1 component of Mandelbrot set */
)$


/* converts angle in radians in range -Pi to Pi
   to turns */
GiveTurn(a):=
( [a,t],
 
  if a<0 then a:a+2*%pi, /* from range (-Pi,Pi) to range (0,2Pi) */
   t:a/(2*%pi) /* from radians to turns */
)$


/* 
http://num.math.uni-goettingen.de/~summer/cheritat.pdf
PARABOLIC IMPLOSION. A MINI-COURSE by ARNAUD CHERITAT

Let Czk be the next term in the expansion of fq:
fq(z) = z + Czk + : : :
We defne attracting axes as the k-1 half lines whose union is the solution to
(C*z^6)/z < 0 ; 
and the repelling axes are for
(C*z^6)/z > 0 ;

*/


GiveTermC(p,q):=
([fq, c, C],
 fq:fn(q,z,c),
 C:coeff(rat(fq),z,q+1),
 c:GiveC(p/q,1), /* parabolic point : rational angle and radius=1.0 */
 C:rectform(float(ev(C)))

)$



GiveTurnOffset(p,q):=
( [C,a,t,offset],
  C:GiveTermC(p,q),
  a:carg(C),
  t:GiveTurn(a),
  tt:t/q, /* turn offset */
  tt:float(rectform(tt))
)$


compile(all);

/* ---------  const -----------  */

p:1;
q:5;
 /* m:exp(2*%i*%pi*p/q); multiplier */


/* --------------- main ----------------------------- */

TurnOffset : GiveTurnOffset(p,q);


/*




// this turn offset is computed with Maxima CAS filename
// star is slightly turned so this offset 
// theory : Complex Dynamics by  Lennart Carleson,Theodore Williams Gamelin  page 40
// TurnOffset = carg(a)/q
// where a = coeff( z^(q+1))

*/

C src code

/*

Adam Majewski
fraktal.republika.pl

c console progam using 
* symmetry
* openMP

 draw  julia sets



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


iHeight =1000
iterMax  = 1000
AR = PixelWidth*15.0; // if you will increase iHeight 
 then increase multiplier or the time of computations will be very long 
real	3m31.345s


File new2001.pgm saved. 
Cx  = 0.356763 
Cy  = 0.328582 
alfax  = 0.154508 
alfay  = 0.475528 
distorsion of image  = 1.000000 

real	407m35.605s


*/



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


// 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 number !!!!!! = (iyMax -iyMin + 1) = iyAboveAxisLength + iyBelowAxisLength +1
// The size of array has to be a positive constant integer 
unsigned int iSize ; // = iWidth*iHeight; 


// memmory 1D array 
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 PixelWidth2; // =  PixelWidth*PixelWidth;
  double PixelHeight; // =(ZyMax-ZyMin)/iYmax;
  double distanceMax;
  double ratio ;
  double TwoPi=2*M_PI;


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

unsigned long int iterMax  = 1000; //iHeight*100;
double ER = 2.0; // Escape Radius for bailout test 
double ER2;
double AR ; // radius around attractor
double AR2; // =AR*AR;


unsigned int period = 5; // period of secondary component joined by root point
  unsigned int denominator;
  double InternalAngle;

unsigned char colors[]={255,230,180, 160,140}; // number of colors >= period


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

  */
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int period)
{
  //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 ( 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 period  there are 2^(period-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 // 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)) 
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;
}




int setup()
{

  
  

  denominator = period;
  InternalAngle = 1.0/((double) denominator);

 

  /* 2D array ranges */
  if (!(iHeight % 2)) iHeight+=1; // it sholud be even number (variable % 2) or (variable & 1)
  iWidth = iHeight;
  iSize = iWidth*iHeight; // size = number 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 ...
  distanceMax = PixelWidth; 
  

  // for numerical optimisation 
  ER2 = ER * ER;
  AR = PixelWidth*15.0; // !!!! important value  ( and precision and time of creation of the pgm image ) 
  AR2= AR*AR;
  PixelWidth2 =  PixelWidth*PixelWidth;
 
  /* 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");

  
 
  return 0;

}



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

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


// forward iteration of complex quadratic polynomial
// fc(z) = z*z +c
// initial point is a z0 = critical point 
// checks if points escapes to infinity
// uses global vars : ER2, Cx, Cy, iterMax
long int GiveLastIteration(double Zx, double Zy )
{ 
    long int iter;
    // z= Zx + Zy*i
    double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
 
     /* initial value of orbit = critical point Z= 0 */
                       
                        Zx2=Zx*Zx;
                        Zy2=Zy*Zy;
                        // for iter from 0 to iterMax because of ++ after last loop
                        for (iter=0; iter<iterMax && ((Zx2+Zy2)<ER2); ++iter)
                        {
                            Zy=2*Zx*Zy + Cy;
                            Zx=Zx2-Zy2 + Cx;
                            Zx2=Zx*Zx;
                            Zy2=Zy*Zy;
                        };
  
   return iter;
}




double GiveTurn(double complex z)
{
 double argument;

 argument = carg(z); //   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
}



/*  */
int GiveNumber(double _Zx0, double _Zy0,double C_x, double C_y, int iMax, double _AR2, double _ZAx, double _ZAy )
{ 
  int i;
  double Zx, Zy; /* z = zx+zy*i */
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  double d, dX, dY; /* distance from z to Alpha  */
  double angle;
 
  Zx=_Zx0; /* initial value of orbit  */
  Zy=_Zy0;
  
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  dX=Zx-_ZAx;
  dY=Zy-_ZAy;
  d=dX*dX+dY*dY;
 
  while ( d>_AR2)
    {
      for(i=0;i<period ;++i) // iMax = period !!!!
       { 
        Zy=2*Zx*Zy + C_y;
        Zx=Zx2-Zy2 +C_x;
        Zx2=Zx*Zx;
        Zy2=Zy*Zy;
       }



      dX=Zx-_ZAx;
      dY=Zy-_ZAy;
      d=dX*dX+dY*dY;
    };

 // divide interior into 5 components 
 angle =GiveTurn(dX +dY*I) - 0.17298227404701; // this turn offset is computed with Maxima CAS filename
// star is slightly turned so this offset 
// theory : Complex Dynamics by  Lennart Carleson,Theodore Williams Gamelin  page 40
// TurnOffset = carg(a)/q
// where a = coeff( z^(q+1))
 if ((0.0<angle) && (angle<0.2)) return 0;
 if ((0.2<angle) && (angle<0.4)) return 1;
 if ((0.4<angle) && (angle<0.6)) return 2;
 if ((0.6<angle) && (angle<0.8)) return 3;
 return 4;
}


 



unsigned char GiveColor(unsigned int ix, unsigned int iy)
{ 
   double Zx, Zy; //  Z= Zx+ZY*i;
  int iter;
  unsigned char color; // gray from 0 to 255 
// unsigned char ColorList[]={255,230,180}; /* shades of gray used in imaage
  //color=0;
  // from screen to world coordinate 
  Zx = GiveZx(ix);
  Zy = GiveZy(iy);
  
 iter  =GiveLastIteration(Zx, Zy );
   if (iter< iterMax ) // point escapes   
   { color = 245; } // exterior 
   else // Filled Julia Set
      { int number =GiveNumber(Zx,Zy, Cx, Cy, period, AR2, creal(alfa), cimag(alfa));
          color = colors[number]; } 


  return color;
}


/* -----------  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) {//printf(" %d from %d\n", ix, ixMax); //info  
                              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;
}



// save data array to pgm file 
int SaveArray2PGMFile( unsigned char data[], unsigned long int iMax)
{
  
  FILE * fp;
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name [10]; /* name of file */
  sprintf(name,"new%lu", iMax); /*  */
  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("Cx  = %f \n", Cx); 
  printf("Cy  = %f \n", Cy);
  // 
 
  printf("alfax  = %f \n", creal(alfa));
  printf("alfay  = %f \n", cimag(alfa));
  printf("iHeght  = %d \n", iHeight);
  
  printf("distorsion of image  = %f \n", ratio);
return 0;
}


/* -----------------------------------------  main   -------------------------------------------------------------*/
int main()
{
// here are procedures for creating image file

  setup();
  c = GiveC(InternalAngle, 1.0, 1) ;
  Cx=creal(c);
  Cy=cimag(c);
  alfa = GiveAlfaFixedPoint(c);
  // compute colors of pixels = image
  //FillArray( data ); // no symmetry
  FillArraySymmetric(data); 
  AddBoundaries(data);
  //  CheckOrientation( data );
  SaveArray2PGMFile(edge , iHeight); // save array (image) to pgm file 
  free(data);
  free(edge);
  info();
  return 0;
}

New c code

/*

Adam Majewski
fraktal.republika.pl

c console progam using 
* symmetry
* openMP

 draw  julia sets



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


iHeight =1000
iterMax  = 1000
AR = PixelWidth*15.0; // if you will increase iHeight 
 then increase multiplier or the time of computations will be very long 
real	3m31.345s


File new2001.pgm saved. 
Cx  = 0.356763 
Cy  = 0.328582 
alfax  = 0.154508 
alfay  = 0.475528 
distorsion of image  = 1.000000 

real	407m35.605s


*/



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


// 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 number !!!!!! = (iyMax -iyMin + 1) = iyAboveAxisLength + iyBelowAxisLength +1
// The size of array has to be a positive constant integer 
unsigned int iSize ; // = iWidth*iHeight; 


// memmory 1D array 
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 PixelWidth2; // =  PixelWidth*PixelWidth;
  double PixelHeight; // =(ZyMax-ZyMin)/iYmax;
  double distanceMax;
  double ratio ;
  double TwoPi=2*M_PI;


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

unsigned long int iterMax  = 1000; //iHeight*100;
double ER = 2.0; // Escape Radius for bailout test 
double ER2;
double AR ; // radius around attractor
double AR2; // =AR*AR;


unsigned int period = 5; // period of secondary component joined by root point
  unsigned int denominator;
  double InternalAngle;

unsigned char colors[]={255,230,180, 160,140}; // number of colors >= period


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

  */
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int period)
{
  //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 ( 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 period  there are 2^(period-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 // 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)) 
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;
}




int setup()
{

  
  

  denominator = period;
  InternalAngle = 1.0/((double) denominator);

 

  /* 2D array ranges */
  if (!(iHeight % 2)) iHeight+=1; // it sholud be even number (variable % 2) or (variable & 1)
  iWidth = iHeight;
  iSize = iWidth*iHeight; // size = number 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 ...
  distanceMax = PixelWidth; 
  

  // for numerical optimisation 
  ER2 = ER * ER;
  AR = PixelWidth*50.0; // !!!! important value  ( and precision and time of creation of the pgm image ) 
  AR2= AR*AR;
  PixelWidth2 =  PixelWidth*PixelWidth;
 
  /* 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");

  
 
  return 0;

}



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

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


// forward iteration of complex quadratic polynomial
// fc(z) = z*z +c
// initial point is a z0 = critical point 
// checks if points escapes to infinity
// uses global vars : ER2, Cx, Cy, iterMax
long int GiveLastIteration(double Zx, double Zy )
{ 
    long int iter;
    // z= Zx + Zy*i
    double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
 
     /* initial value of orbit = critical point Z= 0 */
                       
                        Zx2=Zx*Zx;
                        Zy2=Zy*Zy;
                        // for iter from 0 to iterMax because of ++ after last loop
                        for (iter=0; iter<iterMax && ((Zx2+Zy2)<ER2); ++iter)
                        {
                            Zy=2*Zx*Zy + Cy;
                            Zx=Zx2-Zy2 + Cx;
                            Zx2=Zx*Zx;
                            Zy2=Zy*Zy;
                        };
  
   return iter;
}




double GiveTurn(double complex z)
{
 double argument;

 argument = carg(z); //   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
}



/*  */
int GiveNumber(double _Zx0, double _Zy0,double C_x, double C_y, int iMax, double _AR2, double _ZAx, double _ZAy )
{ 
  int i;
  double Zx, Zy; /* z = zx+zy*i */
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  double d, dX, dY; /* distance from z to Alpha  */
  double angle;
 
  Zx=_Zx0; /* initial value of orbit  */
  Zy=_Zy0;
  
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  dX=Zx-_ZAx;
  dY=Zy-_ZAy;
  d=dX*dX+dY*dY;
 
  while ( d>_AR2)
    {
      for(i=0;i<period ;++i) // iMax = period !!!!
       { 
        Zy=2*Zx*Zy + C_y;
        Zx=Zx2-Zy2 +C_x;
        Zx2=Zx*Zx;
        Zy2=Zy*Zy;
       }



      dX=Zx-_ZAx;
      dY=Zy-_ZAy;
      d=dX*dX+dY*dY;
    };

 // divide interior into 5 components 
 angle =GiveTurn(dX +dY*I); //  - 0.17298227404701; // this turn offset is computed with Maxima CAS filename
// star is slightly turned so this offset 
// theory : Complex Dynamics by  Lennart Carleson,Theodore Williams Gamelin  page 40
// TurnOffset = carg(a)/q
// where a = coeff( z^(q+1))
 if  (angle < 0.16 || angle > 0.97) return 0; // 0.0  < angle < 0.19
 if  (angle < 0.36 ) return 1; // 0.19 < angle
 if  (angle < 0.57 ) return 2; // ((0.4<angle) &&
 if  (angle < 0.775 ) return 3; // ((0.6<angle) && 
 
 return 4; 
}


 



unsigned char GiveColor(unsigned int ix, unsigned int iy)
{ 
   double Zx, Zy; //  Z= Zx+ZY*i;
  int iter;
  unsigned char color; // gray from 0 to 255 
// unsigned char ColorList[]={255,230,180}; /* shades of gray used in imaage
  //color=0;
  // from screen to world coordinate 
  Zx = GiveZx(ix);
  Zy = GiveZy(iy);
  
 iter  =GiveLastIteration(Zx, Zy );
   if (iter< iterMax ) // point escapes   
   { color = 245; } // exterior 
   else // Filled Julia Set
      { int number =GiveNumber(Zx,Zy, Cx, Cy, period, AR2, creal(alfa), cimag(alfa));
          color = colors[number]; } 


  return color;
}


/* -----------  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) {//printf(" %d from %d\n", ix, ixMax); //info  
                              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;
}



// save data array to pgm file 
int SaveArray2PGMFile( unsigned char data[], unsigned long int iMax)
{
  
  FILE * fp;
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name [10]; /* name of file */
  sprintf(name,"new%lu", iMax); /*  */
  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("Cx  = %f \n", Cx); 
  printf("Cy  = %f \n", Cy);
  // 
 
  printf("alfax  = %f \n", creal(alfa));
  printf("alfay  = %f \n", cimag(alfa));
  printf("iHeght  = %d \n", iHeight);
  
  printf("distorsion of image  = %f \n", ratio);
return 0;
}


/* -----------------------------------------  main   -------------------------------------------------------------*/
int main()
{
// here are procedures for creating image file

  setup();
  c = GiveC(InternalAngle, 1.0, 1) ;
  Cx=creal(c);
  Cy=cimag(c);
  alfa = GiveAlfaFixedPoint(c);
  // compute colors of pixels = image
  //FillArray( data ); // no symmetry
  FillArraySymmetric(data); 
  AddBoundaries(data);
  //  CheckOrientation( data );
  SaveArray2PGMFile(edge , iHeight); // save array (image) to pgm file 
  free(data);
  free(edge);
  info();
  return 0;
}

result of new code :

memory is OK 
axis of symmetry 
 find boundaries in data array using  Sobel filter
File new2001.pgm saved. 
Cx  = 0.356763 
Cy  = 0.328582 
alfax  = 0.154508 
alfay  = 0.475528 
iHeght  = 2001 
distorsion of image  = 1.000000 

real	0m7.806s
user	1m1.933s
sys	0m0.042s

Rererences

File history

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

	Date/Time	Thumbnail	Dimensions	User	Comment
current	16:07, 5 December 2025		1,993 × 1,993 (33 KB)	wikimediacommons>Sebastian Wallroth	Cropped < 1 % horizontally, < 1 % vertically, 1 % areawise using CropTool with precise mode.

File usage

The following 2 pages use this file:

File:Parabolic Julia set for internal angle 1 over 5 - new method.png

Summary

Algorithm

C src code

Rererences

Captions

Items portrayed in this file

depicts

creator

some value

copyright status

copyrighted

copyright license

Creative Commons Attribution-ShareAlike 3.0 Unported

GNU Free Documentation License, version 1.2 or later

inception

26 December 2012

source of file

original creation by uploader

File history

File usage

Navigation menu

File:Parabolic Julia set for internal angle 1 over 5 - new method.png

Summary

Algorithm

C src code

Rererences

Captions

Items portrayed in this file

some value

26 December 2012

File history

File usage

Navigation menu

Search