Fractal Related Discussion > Fractal Mathematics And New Theories

 rational function

(1/16) > >>

Adam Majewski:
I want to draw   Julia set for rational function (image bo36.png)  by  Michael Becker : 

--- Quote ---f(z)=1/((0,15+0,15i)z5+z3+(-3+3i)z), dargestellt auf [-3;3]x[-3;3].
--- End quote ---

It look like connected Julia set. To speed up computaions I want to find critical points and periodic points


--- Code: ---
kill(all);
remvalue(all);
display2d:false;


f(z):= 1/(a*z^5+ z^3 + b*z);


/*
converts complex number z = x*y*%i
to the list in a draw format: 
[x,y]
*/
d(z):=[float(realpart(z)), float(imagpart(z))]$

ToPoints(myList):= points(map('d,myList))$


/* give Draw List from one point*/
ToPoint(z):=points([d(z)])$





/* f(z) is used as a global function
   I do not know how to put it as a argument */

GiveOrbit(z0):=
block(
 [z,Orbit, iLength ],
 iLength : 100,
 z:z0,
 Orbit:[z0],
 for i:1 thru iLength step 1 do
        ( z:f(z),
          z: float(z),
          z:rectform(z),
          Orbit:endcons(z,Orbit)),
 
 Orbit:ToPoints(Orbit),       
 return(Orbit)

)$




/* find fixed points  returns a list */
GiveFixedPoints():= block
(
  [s,z],
  s:solve(f(z)=z),
  /* remove "z="  from list s */
  s:map('rhs,s),
  s:map('rectform,s),
  s:map('float,s),
  return(s)
)$




/*compile(all);*/


a: (0.15+0.15*%i);
b: (-3+3*%i);


dz: diff(f(z) ,z,1)$

dz:ratsimp(float(ev(dz)));
 
 /*
 dz = -((150*%i+150)*z^4+600*z^2+600*%i-600) /(9*%i*z^10+(60*%i+60)*z^8-160*z^6+(1200*%i-1200)*z^4-3600*%i*z^2)
*/

s: solve(dz=0)$
s:map('rhs,s)$
s:map('float,s)$
s:map('rectform,s)$
s:map('float,s);
/*


 [0.389374737779177*%i-0.9400337727919567,
        0.9400337727919567-0.389374737779177*%i,
        (-1.816005797447402*%i)-0.7522142306508819,
        1.816005797447402*%i+0.7522142306508819]
 
*/


/* r1:solve(f(z)=z); */




iLength:100;
orbit1 :GiveOrbit(s[1])$
orbit2 :GiveOrbit(s[2])$
orbit3 :GiveOrbit(s[3])$
orbit4 :GiveOrbit(s[4])$



path:"~/Dokumenty/ijon/3/36/"$  /*  if empty then file is in a home dir */

load(draw); /* ( interface to gnuplot ) by Mario Rodriguez Riotorto http://www.telefonica.net/web2/biomates */

draw2d(
    title = concat("Critical orbit for f(z)=z^2 +", string(c)),
    terminal  = png,
    user_preamble = "set size square", /*    */
    file_name = concat(path ,"8"),
    dimensions = [1500,1500],    /* Since Maxima 5.23, pic_width and pic_height are deprecated. */
    xrange = [-3,3],
    yrange = [-3,3],
    xlabel     = "z.re ",
    ylabel     = "z.im",

   
   
    point_type    = 7,
    points_joined = false,
    point_size    = 0.8,
    key=" critical orbit ",
    color             =red,
    orbit1,
    orbit2,
    orbit3,
    orbit4,
   
    point_size    = 1.2,
    key="",
    ToPoint(s[1]),
    ToPoint(s[2]),
    ToPoint(s[3]),
    key= "critical point",
    color           = blue,
    ToPoint(s[4])
   
    /*
   

    key= "fixed point",
    color           = black,
    alfa
*/
 );


--- End code ---

My image of critical orbits  shows that critical orbits tend to z=0, but on the image it seems to be in the exterior
How to find periodic points ?







Linkback: https://fractalforums.org/index.php?topic=4279.0

claude:

--- Quote from: Adam Majewski on June 17, 2021, 05:16:22 PM ---critical orbits tend to z=0

--- End quote ---
Are you sure about that?  Possibly tending to fixed points near (but not equal to) 0?  I have not investigated, but it looks like there is a gap in the middle, in both images.

EDIT fixed points and period 2 points:

--- Code: ---(%i10) ToPoints(solve(z = f(z), z));

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15
(%o10) points([[z = - 0.1728968618954043, 0.0 = 0.4517976427751508],
[z = 0.1728968618954043, 0.0 = - 0.4517976427751508],
[z = - 1.542012716025067, 0.0 = 0.5826223205347855],
[z = 1.542012716025067, 0.0 = - 0.5826223205347855],
[z = - 1.052616733765589, 0.0 = - 2.511067897204381],
[z = 1.052616733765589, 0.0 = 2.511067897204381]])
(%i18) map('allroots, solve(f(f(z))=z));

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15

rat: replaced 0.15 by 3/20 = 0.15
(%o18) [[z = 0.0], [z = 0.1004507524894826 - 0.04160806403172667 %i,
z = 0.04160806403172667 %i - 0.1004507524894826,
z = (- 0.1972126334382917 %i) - 0.4417257204225227,
z = 0.1728968618954057 - 0.4517976427751519 %i,
z = 0.1972126334382916 %i + 0.4417257204225231,
z = 0.4517976427751519 %i - 0.172896861895406,
z = (- 2.008427151852299 %i) - 0.8319177653355924,
z = 1.486804283840914 - 0.6158544989534238 %i,
z = 0.5826223205446455 %i - 1.542012716044921,
z = 0.6783914544446792 %i - 1.502343841878277,
z = 1.545643025048946 - 0.6402263035701801 %i,
z = 1.542012716021183 - 0.5826223205577208 %i,
z = 0.6158544989143758 %i - 1.48680428382238,
z = 0.6402263036121009 %i - 1.545643025068914,
z = 1.502343841903455 - 0.678391454434491 %i,
z = 2.511067896520314 %i + 1.052616752123668,
z = 2.008427151854269 %i + 0.83191776533551,
z = (- 2.519905568763936 %i) - 1.031280709377238,
z = (- 2.511067897025448 %i) - 1.052616732010429,
z = 2.519905568548988 %i + 1.03128068923336,
z = 2.496853258644555 %i + 1.034230487296836,
z = (- 2.532331263153173 %i) - 1.048925954202523,
z = (- 2.496853270435978 %i) - 1.034230487368577,
z = 2.532331275662722 %i + 1.048925954304978]]

--- End code ---

Adam Majewski:

--- Quote from: claude on June 17, 2021, 05:27:02 PM ---Are you sure about that?  Possibly tending to fixed points near (but not equal to) 0?  I have not investigated, but it looks like there is a gap in the middle, in both images.

--- End quote ---

Ok. It looks that z=0 will be inside ( ?center) of very small component which is smth like limit of white components.

So period 1 point will be not helpfull for me.

I ask about it because usually testing for attracting periodic points makes the code faster.

claude:
see edit, period 2 points include 0 and +/-(0.04160806403172667 %i - 0.1004507524894826)

Adam Majewski:

--- Code: ---
/*

https://web.archive.org/web/20161024194536/http://www.ijon.de/mathe/julia/some_julia_sets_3.html

a: (0,15+0,15i)
b: (-3+3i)
f: 1/(a*z^5+ b*z^3 + z);
 dz = -(5*a*z^4+3*b*z^2+1)/(a*z^5+z^3+b*z)^2


  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero
 
 
 
  Structure of a program or how to analyze the program
 
 
  ============== Image X ========================
 
  DrawImageOf -> DrawPointOf -> ComputeColorOf ( FunctionTypeT FunctionType , complex double z) -> ComputeColor
 
 
  check only last function  which computes color of one pixel for given Function Type
 
 

   
  ==========================================

 
  ---------------------------------
  indent d.c
  default is gnu style
  -------------------



  c console progam
 
  export  OMP_DISPLAY_ENV="TRUE"
  gcc d.c -lm -Wall -march=native -fopenmp
  time ./a.out > b.txt


  gcc d.c -lm -Wall -march=native -fopenmp


  time ./a.out

  time ./a.out >i.txt
  time ./a.out >e.txt
 
 
 
 
 
 
  convert -limit memory 1000mb -limit disk 1gb dd30010000_20_3_0.90.pgm -resize 2000x2000 10.png

 
 
 
*/

#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
#include <limits.h> // Maximum value for an unsigned long long int



// https://sourceforge.net/p/predef/wiki/Standards/

#if defined(__STDC__)
#define PREDEF_STANDARD_C_1989
#if defined(__STDC_VERSION__)
#if (__STDC_VERSION__ >= 199409L)
#define PREDEF_STANDARD_C_1994
#endif
#if (__STDC_VERSION__ >= 199901L)
#define PREDEF_STANDARD_C_1999
#endif
#endif
#endif




/* --------------------------------- global variables and consts ------------------------------------------------------------ */


//FunctionType
typedef enum  {Fatou = 0, IntLSM =1 , ExtLSM = 2, LSM = 3, DEM = 4, Unknown = 5 , BD = 6, MBD = 7 , SAC = 8, DLD = 9, ND = 10 , NP= 11, POT = 12 , Blend = 13

} FunctionTypeT;
// FunctionTypeT FunctionType;

// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1
//unsigned int ix, iy; // var
static unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
static unsigned int ixMax; //
static unsigned int iWidth; // horizontal dimension of array

static unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iyMax; //

static unsigned int iHeight = 1000; // 
// The size of array has to be a positive constant integer
static unsigned long long int iSize; // = iWidth*iHeight;

// memmory 1D array
unsigned char *data;
unsigned char *edge;
//unsigned char *edge2;

// unsigned int i; // var = index of 1D array
//static unsigned int iMin = 0; // Indexes of array starts from 0 not 1
static 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



// see SetPlane

double radius = 3.0;
complex double center = 0.0 ;
double  DisplayAspectRatio  = 1.0; // https://en.wikipedia.org/wiki/Aspect_ratio_(image)
// dx = dy compare setup : iWidth = iHeight;
double ZxMin; //= -1.3; //-0.05;
double ZxMax;// = 1.3; //0.75;
double ZyMin;// = -1.3; //-0.1;
double ZyMax;// = 1.3; //0.7;
double PixelWidth; // =(ZxMax-ZxMin)/ixMax;
double PixelHeight; // =(ZyMax-ZyMin)/iyMax;

// dem
double BoundaryWidth ; //= 1.0*iWidth/2000.0  ; //  measured in pixels ( when iWidth = 2000)
double distanceMax ; //= BoundaryWidth*PixelWidth;


double ratio;


/*
  ER = pow(10,ERe);
  AR = pow(10,-ARe);
*/
//int ARe ; // increase ARe until black ( unknown) points disapear
//int ERe ;
double ER;
double ER2; //= 1e60;
double AR; // bigger values do not works
double AR2;

double AR_max;
//double AR12;



int IterMax = 100000;
int IterMax_LSM = 1000;
int IterMax_DEM = 100000;

/* colors = shades of gray from 0 to 255

   unsigned char colorArray[2][2]={{255,231},    {123,99}};
   color = 245;  exterior
*/
unsigned char iColorOfExterior = 245;
unsigned char iColorOfInterior = 99;
//unsigned char iColorOfInterior2 = 183;
unsigned char iColorOfBoundary = 0;
unsigned char iColorOfUnknown = 5;

// pixel counters
unsigned long long int uUnknown = 0;
unsigned long long int uInterior = 0;
unsigned long long int uExterior = 0;



// critical points




// critical points
//complex double zcr = -0.35*I; // only one critical point
//complex double zc2 = -2.2351741790771484375e-08+9.4296410679817199707e-09*I;


/*

  0.389374737779177*%i-0.9400337727919567,
        0.9400337727919567-0.389374737779177*%i,
        (-1.816005797447402*%i)-0.7522142306508819,
        1.816005797447402*%i+0.7522142306508819]

*/

const complex double zp2 = 0.389374737779177*I-0.9400337727919567; //-0.04160806403172667*I - 0.1004507524894826;

/*
a: (0.15+0.15*%i);
b: (-3+3*%i);
f: 1/(a*z^5+ b*z^3 + z);
 dz = -(5*a*z^4+3*b*z^2+1)/(a*z^5+b*z^3+z)^2

*/
const complex double a =  0.15 + 0.15*I;
const complex double b = -3.0  + 3.0*I;
//const int period = 1;

// periodic points = attractors
//complex double z1 =  0.0 ; //fixed point (period  1) =  attracting cycle



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


// complex function
complex double f(const complex double z0) {

double complex z = z0;
z = 1/(a*z*z*z*z*z + z*z*z + b*z);
return  z;
}


//
complex double dfz(const complex double z0) {


// dz= -(5*a*z^4+3*z^2+b) / (a*z^5+z^3+b*z)^2
double complex z = z0;
complex double z2= z*z;
complex double numerator = -5.0*a*z2*z2 - 3.0*z2 - b ;
complex double denom = a*z2*z2*z + z2*z + b*z;
denom = denom*denom;


return  numerator/denom;

}









//------------------complex numbers -----------------------------------------------------





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

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


complex double
GiveZ (int ix, int iy)
{
  double Zx = GiveZx (ix);
  double Zy = GiveZy (iy);

  return Zx + Zy * I;




}



double cabs2(complex double z){

  return creal(z)*creal(z)+cimag(z)*cimag(z);


}



// ****************** DYNAMICS = trap tests ( target sets) ****************************


/* -----------  array functions = drawing -------------- */

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















//
unsigned char ComputeColorOfFatou (complex double z)
{





  double r2;


  int i; // number of iteration
  for (i = 0; i < IterMax; ++i)
    {




        r2 =cabs2(z);

      if (r2 > ER2) // esaping = exterior
{
  uExterior += 1;
  return iColorOfExterior;
}


if (cabs2(zp2-z) < AR2) //
{
 
  return iColorOfInterior;
}


     
z = f(z); // complex iteration f(z)=z^3 + c

    }

 
  return iColorOfUnknown;


}









// f(z)=1+z−3z2−3.75z3+1.5z4+2.25z5
unsigned char ComputeColorOfLSM (complex double z)
{





  double r2;


  int i; // number of iteration
  for (i = 0; i < IterMax_LSM; ++i)
    {




      // complex iteration f(z)=z^3 + c
            r2 =cabs2(z);

      if (r2 > ER2) // esaping = exterior
{
  uExterior += 1;
  return 255- ((i*15) % 255);
}

       if (cabs2(zp2-z) < AR2)  // interior
{
 
  return iColorOfInterior;
}

z = f(z);

    }

  return iColorOfUnknown;


}













// ***************************************************************************************************************************
// ************************** DEM/J*****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfDEMJ(complex double z){
  // https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#DEM.2FJ


 
  int nMax = IterMax_DEM;
  complex double dz = 1.0; //  is first derivative with respect to z.
  double distance;
  double cabsz;

  int n;

  for (n=0; n < nMax; n++){ //forward iteration
    cabsz = cabs(z);
    if (cabsz > 1e60 || cabs(dz)> 1e60) break; // big values
    if (cabs2(zp2-z) < AR2)  return iColorOfInterior; // falls into finite attractor = interior
 
    dz = dfz(z)*dz;
    z = f(z) ; /* forward iteration : complex cubic polynomial */
  }
 
 
  distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
  if (distance <distanceMax) return iColorOfBoundary; // distanceMax = BoundaryWidth*PixelWidth;
  // else
 
  return iColorOfExterior;

 
}






/* ==================================================================================================
   ============================= Draw functions ===============================================================
   =====================================================================================================
*/
unsigned char ComputeColor(FunctionTypeT FunctionType, complex double z){

  unsigned char iColor;



  switch(FunctionType){
 
  case Fatou :{iColor = ComputeColorOfFatou(z); break;}
 
 
  // case IntLSM :{iColor = ComputeColorOfIntLSM(z); break;}

  // case ExtLSM :{iColor = ComputeColorOfExtLSM(z); break;}
 
  case LSM :{iColor = ComputeColorOfLSM(z); break;}
 
 

  case DEM : {iColor = ComputeColorOfDEMJ(z); break;}

/*
  case Unknown : {iColor = ComputeColorOfUnknown(z); break;}

  case BD : {iColor = ComputeColorOfBD(z); break;}

  case MBD : {iColor = ComputeColorOfMBD(z); break;}

  case SAC : {iColor = ComputeColorOfSAC(z); break;}
 
  case DLD : {iColor = ComputeColorOfDLD(z); break;}

  case ND : {iColor = ComputeColorOfND(z); break;}

  case NP : {iColor = ComputeColorOfNP(z); break;}

  case POT : {iColor = ComputeColorOfPOT(z); break;}

  case Blend : {iColor = ComputeColorOfBlend(z); break;}
  */

  default: {}


  }

  return iColor;



}


// plots raster point (ix,iy)
int DrawPoint ( unsigned char A[], FunctionTypeT FunctionType, int ix, int iy)
{
  int i; /* index of 1D array */
  unsigned char iColor;
  complex double z;


  i = Give_i (ix, iy); /* compute index of 1D array from indices of 2D array */
 
  z = GiveZ(ix,iy);
 

  iColor = ComputeColor(FunctionType, z);
  A[i] = iColor ; //
 
  return 0;
}




int DrawImage ( unsigned char A[], FunctionTypeT FunctionType)
{
  unsigned int ix, iy; // pixel coordinate

  fprintf (stdout, "compute image \n");
  // for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax, uUnknown, uInterior, uExterior)
  for (iy = iyMin; iy <= iyMax; ++iy)
    {
      fprintf (stdout, " %d from %d \r", iy, iyMax); //info
      for (ix = ixMin; ix <= ixMax; ++ix)
DrawPoint(A, FunctionType, ix, iy); // 
    }
  fprintf (stdout, "\n"); //info
  return 0;
}







int PlotPoint(complex double z, unsigned char A[]){


  unsigned int ix = (creal(z)-ZxMin)/PixelWidth;
  unsigned int iy = (ZyMax - cimag(z))/PixelHeight;
  unsigned int i = Give_i(ix,iy); /* index of _data array */


  A[i]= 0; //255-A[i]; // Mark point with inveres color


  return 0;

}






// ***********************************************************************************************
// ********************** edge detection usung Sobel filter ***************************************
// ***************************************************************************************************

// from Source to Destination
int ComputeBoundaries(unsigned char S[], unsigned char D[])
{
 
  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;
  // boundaries are in D  array ( global var )
 
  // clear D array
  memset(D, iColorOfExterior, iSize*sizeof(*D)); // for heap-allocated arrays, where N is the number of elements = FillArrayWithColor(D , iColorOfExterior);
 
  // printf(" find boundaries in S array using  Sobel filter\n");   
#pragma omp parallel for schedule(dynamic) private(i,iY,iX,Gv,Gh,G) shared(iyMax,ixMax)
  for(iY=1;iY<iyMax-1;++iY){
    for(iX=1;iX<ixMax-1;++iX){
      Gv= S[Give_i(iX-1,iY+1)] + 2*S[Give_i(iX,iY+1)] + S[Give_i(iX-1,iY+1)] - S[Give_i(iX-1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[Give_i(iX+1,iY-1)];
      Gh= S[Give_i(iX+1,iY+1)] + 2*S[Give_i(iX+1,iY)] + S[Give_i(iX-1,iY-1)] - S[Give_i(iX+1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[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) {D[i]=255;} /* background */
      else {D[i]=0;}  /* boundary */
    }
  }
 
   
 
  return 0;
}



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

















// *******************************************************************************************
// ********************************** save A array to pgm file ****************************
// *********************************************************************************************

int SaveArray2PGMFile (unsigned char A[],  char * n, char *comment)
{

  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 */
  snprintf (name, sizeof name, "%s", n ); /*  */
  char *filename = strcat (name, ".pgm");
  char long_comment[200];
  sprintf (long_comment, "fc(z)=z^2+z*c  %s", comment);





  // save image array 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", long_comment, iWidth, iHeight, MaxColorComponentValue); // write header to the file
  size_t rSize = fwrite (A, sizeof(A[0]), iSize, fp); // write whole array with image data bytes to the file in one step
  fclose (fp);

  // info
  if ( rSize == iSize)
  {
  printf ("File %s saved ", filename);
  if (long_comment == NULL || strlen (long_comment) == 0)
    printf ("\n");
  else { printf (". Comment = %s \n", long_comment); }
  }
  else {printf("wrote %zu elements out of %llu requested\n", rSize,  iSize);}

  return 0;
}




int PrintCInfo ()
{

  printf ("gcc version: %d.%d.%d\n", __GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__); // https://stackoverflow.com/questions/20389193/how-do-i-check-my-gcc-c-compiler-version-for-my-eclipse
  // OpenMP version is displayed in the console : export  OMP_DISPLAY_ENV="TRUE"

  printf ("__STDC__ = %d\n", __STDC__);
  printf ("__STDC_VERSION__ = %ld\n", __STDC_VERSION__);
  printf ("c dialect = ");
  switch (__STDC_VERSION__)
    { // the format YYYYMM
    case 199409L:
      printf ("C94\n");
      break;
    case 199901L:
      printf ("C99\n");
      break;
    case 201112L:
      printf ("C11\n");
      break;
    case 201710L:
      printf ("C18\n");
      break;
      //default : /* Optional */

    }

  return 0;
}


int
PrintProgramInfo ()
{


  // display info messages
  printf ("Numerical approximation of Julia set for F(z) =   z*c) \n");
  //printf ("parameter C = ( %.16f ; %.16f ) \n", creal (C), cimag (C));
 

  printf ("Image Width = %f in world coordinate\n", ZxMax - ZxMin);
  printf ("PixelWidth = %.16f \n", PixelWidth);
  //printf ("AR = %.16f = %f *PixelWidth = %f %% of ImageWidth \n", AR, AR / PixelWidth, AR / ZxMax - ZxMin);



  // image corners in world coordinate
  // center and radius
  // center and zoom
  // GradientRepetition
  printf ("Maximal number of iterations = iterMax = %d \n", IterMax);
  printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  //




  return 0;
}



int SetPlane(complex double center, double radius, double a_ratio){

  ZxMin = creal(center) - radius*a_ratio;
  ZxMax = creal(center) + radius*a_ratio; //0.75;
  ZyMin = cimag(center) - radius; // inv
  ZyMax = cimag(center) + radius; //0.7;
  return 0;

}



// Check Orientation of z-plane image : mark first quadrant of complex plane
// it should be in the upper right position
// uses global var :  ...
int CheckZPlaneOrientation(unsigned char A[] )
{
 
double Zx, Zy; //  Z= Zx+ZY*i;
unsigned i; /* index of 1D array */
unsigned int ix, iy; // pixel coordinate

fprintf(stderr, "compute image CheckOrientation\n");
  // for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy, i, Zx, Zy) shared(A, ixMax , iyMax)
for (iy = iyMin; iy <= iyMax; ++iy){
    //fprintf (stderr, " %d from %d \r", iy, iyMax); //info
    for (ix = ixMin; ix <= ixMax; ++ix){
    // from screen to world coordinate
    Zy = GiveZy(iy);
    Zx = GiveZx(ix);
  i = Give_i(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;
}







// *****************************************************************************
//;;;;;;;;;;;;;;;;;;;;;;  setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
// **************************************************************************************

int setup ()
{

  fprintf (stderr, "setup start\n");






  /* 2D array ranges */

  iWidth = iHeight* DisplayAspectRatio ;
  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].
  //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].

 
  SetPlane( center, radius,  DisplayAspectRatio );
  /* Pixel sizes */
  PixelWidth = (ZxMax - ZxMin) / ixMax; //  ixMax = (iWidth-1)  step between pixels in world coordinate
  PixelHeight = (ZyMax - ZyMin) / iyMax;
  ratio = ((ZxMax - ZxMin) / (ZyMax - ZyMin)) / ((double) iWidth / (double) iHeight); // it should be 1.000 ...

  ER = 200.0; //
  ER2 = ER*ER;
  //AR_max = 5*PixelWidth*iWidth/2000.0 ; // adjust first number
  // GiveTunedAR(const int i_Max, const complex double zcr, const double c, const double zp){
  AR = 0.003; // GiveTunedAR(200, C, AR_max);
  AR2 = AR * AR;
  //AR12 = AR/2.0;
 
 
    BoundaryWidth = 1.0*iWidth/2000.0  ; //  measured in pixels ( when iWidth = 2000)
  distanceMax = BoundaryWidth*PixelWidth;



  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc (iSize * sizeof (unsigned char));

  edge = malloc (iSize * sizeof (unsigned char));
  if (data == NULL || edge == NULL)
    {
      fprintf (stderr, " Could not allocate memory");
      return 1;
    }





 


  fprintf (stderr, " end of setup \n");

  return 0;

} // ;;;;;;;;;;;;;;;;;;;;;;;;; end of the setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




int end ()
{


  fprintf (stderr, " allways free memory (deallocate )  to avoid memory leaks \n"); // https://en.wikipedia.org/wiki/C_dynamic_memory_allocation
  free (data);
  free(edge);


  PrintProgramInfo ();
  PrintCInfo ();
  return 0;

}

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

int main ()
{
  setup ();

   DrawImage (data, Fatou); // first
  SaveArray2PGMFile (data,  "Fatou_1" , "Fatou ");
 
  DrawImage (data, LSM); // first
  SaveArray2PGMFile (data,  "LSM_1" , "LSM ");
 
 
 
  DrawImage (data, DEM); // first
  SaveArray2PGMFile (data,  "DEM_1" , "DEM ");
 

 
  end ();

  return 0;
}




--- End code ---

Navigation

[0] Message Index

[#] Next page

Go to full version
Powered by SMFPacks Rates (Facepunch) Mod
Powered by SMFPacks SEO Pro Mod | Sitemap
Powered by Advanced Topic Prefix Pro
Powered by SMFPacks Advanced Attachments Uploader Mod