File:Julia set p(z)= z^3+(1.0149042485835864102+0.10183008497976470119i)*z; (zoom).png
Jump to navigation
Jump to search
Size of this preview: 600 × 600 pixels. Other resolutions: 240 × 240 pixels | 480 × 480 pixels | 768 × 768 pixels | 1,024 × 1,024 pixels | 2,000 × 2,000 pixels.
Original file (2,000 × 2,000 pixels, file size: 1.65 MB, MIME type: image/png)
Captions
Captions
Julia set p(z)= z^3+(1.0149042485835864102+0.10183008497976470119i)*z. Zoom and critical orbit
Summary
[edit]| Description |
English: Julia set p(z)= z^3+(1.0149042485835864102+0.10183008497976470119i)*z. Zoom and critical orbit. Location by Michael Becker[1] Image made with LCM ( level curves for interior and exterior). Level curves cross at critical point and its preimages. One can see spiral from attracting fixed point to repelling fixed point ( z= 0) which is a place with high density of level curves. Points of critical orbit ( including crirital point and attractor) are on the level curves like notes on the musical staff Deutsch: f(z)=z3+dz mit d=1,02*e0,1i, dargestellt auf [-0.5;0.5]x[-0.5;0.5]. The Julia set (boundary of filled-n Julia set) itself is not drawn: we see it as the locus of points where the level curves (= the boundaries of level sets) are especially close to each other = a place with high density of level curves. |
| Date | |
| Source | Own work |
| Author | Adam majewski |
| Other versions |
|
%3D_z%5E3%2B(1.0149042485835864102%2B0.10183008497976470119i)*z.png)
_%3D_z%5E3_%2Bz*(0.1008317508132964*i_%2B_1.004954206930806).png)
_%3D_z%5E2_%2B0.25.svg)
%3D_z%5E3%2B(1.0149042485835864102%2B0.10183008497976470119i)*z;_(zoom).png&action=edit§ion=1){kind=link}
Licensing
[edit]%3D_z%5E3%2B(1.0149042485835864102%2B0.10183008497976470119i)*z;_(zoom).png&action=edit§ion=2){kind=link}
I, the copyright holder of this work, hereby publish it under the following license:
This file is licensed under the Creative Commons Attribution-Share Alike 4.0 International license.
- You are free:
- to share – to copy, distribute and transmit the work
- to remix – to adapt the work
- Under the following conditions:
- attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
- share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.
c source code
[edit]%3D_z%5E3%2B(1.0149042485835864102%2B0.10183008497976470119i)*z;_(zoom).png&action=edit§ion=3){kind=link}
/*
Adam Majewski
adammaj1 aaattt o2 dot pl // o like oxygen not 0 like zero
console program in c programing language
===============================================================
coefficients read from input file ijon_b011.txt
degree 3 coefficient = ( +1.0000000000000000 +0.0000000000000000*i)
degree 1 coefficient = ( +1.0149042485835864 +0.1018300849797647*i)
Input polynomial p(z)=(1+0i)*z^3+(1.0149042485835864102+0.10183008497976470119i)*z^1
derivative dp/dz = (3+0i)*z^2+(1.0149042485835864102+0.10183008497976470119i)
2 critical points found
cp#0: 0.029142613176165576422,-0.58236647233272664792 . It's critical orbit is bounded and enters cycle #0 length=1 and it's stability = |multiplier|=0.99134 =parabolic
internal angle = 0.96706874577169499307
cycle = {
0.20977461555364212975,-0.24271307734496930242 ; }
cp#1: -0.029142613176165569483,0.58236647233272664792 . It's critical orbit is bounded and enters cycle #1 length=1 and it's stability = |multiplier|=0.99134 =parabolic
internal angle = 0.96706874577169499307
cycle = {
-0.20977461555364212975,0.24271307734496930242 ; }
==============================================
Structure of a program or how to analyze the program
============== Image of (FunctionType, PlaneInversion) ========================
DrawImageOf(X,y) -> DrawPoint(X,y) -> ComputeColorOf(X,y)
so check only compute color function
==========================================
---------------------------------
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 >a.txt
./g.sh
make
----------------------
real 0m19,809s
user 2m26,763s
sys 0m0,161s
*/
#include <stdio.h>
#include <stdlib.h> // malloc
#include <string.h> // strcat
#include <math.h> // M_PI; needs -lm also
#include <complex.h> // complex numbers : https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
#include <omp.h> // OpenMP
/* --------------------------------- global variables and consts ------------------------------------------------------------ */
//
double complex c; // parameter of function fc(z)=z^3 + c*z
// attractors
double complex za0 = 0.20977461555364212975 -0.24271307734496930242 *I; // first attracting point of period 1 = cycle #0
double complex za1 = -0.20977461555364212975 +0.24271307734496930242 *I; // second attracting point of period 1 = cycle #1
// for bd
double za0im;
double za1im;
double complex zcr0 = 0.029142613176165576422-0.58236647233272664792 *I; // first critical = cycle #0
double complex zcr1 = -0.029142613176165569483+0.58236647233272664792 *I; // second critical point = cycle #1
int NumberOfImages = 0;
//FunctionType
typedef enum {Fatou, IntLSM, ExtLSM , LSM, DEM, Unknown, BD, MBD , SAC, DLD, ND, NP, POT, Blend
} FunctionTypeT;
// FunctionTypeT FunctionType;
int PlaneInversion = 0; // boolean 1 = w = 1/z plane; 0 = z plane
// 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 = 20000; //
// The size of array has to be a positive constant integer
static unsigned int iSize; // = iWidth*iHeight;
// ----------memmory 1D arrays ==================
// unsigned char = for 1 byte ( 8 bit) colors
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 SetZPlane
double radius = 0.5;
complex double center = 0.0;
double DisplayAspectRatio = 1.0; // https://en.wikipedia.org/wiki/Aspect_ratio_(image)
// z plane = dynamic plane
double ZxMin ; //-0.05;
double ZxMax ; //0.75;
double ZyMin ; //-0.1;
double ZyMax ; //0.7;
double PixelWidth; // =(ZxMax-ZxMin)/ixMax;
double PixelHeight; // =(ZyMax-ZyMin)/iyMax;
double ratio;
// w plane = 1/z plane
double WxMin = -2000; //-0.05;
double WxMax = 2000; //0.75;
double WyMin = -2000; //-0.1;
double WyMax = 2000; //0.7;
double wPixelWidth; // =(WxMax-WxMin)/ixMax;
double wPixelHeight; // =(WyMax-WyMin)/iyMax;
static unsigned long int iterMax = 1000000; //iHeight*100;
const long int iterMax_LSM = 1000;
const int iterMax_DLD = 200; // N in wiki = fixed number : maximal number of iterations
const int iterMax_pot = 400; // potential
int iterMax_DEM = 1000;
double ER = 200.0; // EscapeRadius for bailout test
double EscapeRadius=1000000; // = ER big !!!!
double ER_LSM ; // see GiveER_LSM // 27.764 = manually find value such that level curves of escape time cross critical point and it's preimages
double ER_DLD ; // see GiveER_LSM // 27.764 = manually find value such that level curves of escape time cross critical point and it's preimages
double ER_NP = 100.0;
double ER_pot = 100000.0; // sqrt(1e24)
double ER_DEM = 1000.0; // sqrt(1e24)
double loger; // = log(ER_LSM); // for texture
static double TwoPi=2.0*M_PI; // texture
double MaxFinalRadius;
double AR; // Attracting Radius = radius of circle around attractor ( trap)
double AR_max;
// SAC/J
double lnER; // ln(ER)
int i_skip = 2; // exclude (i_skip+1) elements from average
unsigned char s = 7; // stripe density
double BoundaryWidth = 3.0; // % of image width
double distanceMax; //distanceMax = BoundaryWidth*PixelWidth;
// ------------- DLD ----------------------
double p = 0.180; //0.01444322; //
// DLD colors
//double me = 1.0;
double mi = 0.9;
// potential
double MaxImagePotential;
/* colors = shades of gray from 0 to 255 */
unsigned char iColorOfExterior = 250;
unsigned char iColorOfInterior = 200;
unsigned char iColorOfInterior1 = 210;
unsigned char iColorOfInterior2 = 180;
unsigned char iColorOfBoundary = 0;
unsigned char iColorOfUnknown = 30;
/* ------------------------------------------ functions -------------------------------------------------------------*/
/**
* Find maximum between two numbers.
https://codeforwin.org/2016/02/c-program-to-find-maximum-and-minimum-using-functions.html
*/
double max(double n1, double n2)
{
return (n1 > n2 ) ? n1 : n2;
}
//---------------------
double min(double n1, double n2)
{
return (n1 < n2 ) ? n1 : n2;
}
double clip(double d){
return (d> 1.0) ? 1.0 : d;
}
double frac(double d){
double fraction = d - ((long)d);
return fraction;
}
//------------------complex numbers -----------------------------------------------------
double c_arg(complex double z)
{
double arg;
arg = carg(z);
if (arg<0.0) arg+= TwoPi ;
return arg;
}
double c_turn(complex double z)
{
double arg;
arg = c_arg(z);
return arg/TwoPi;
}
// 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;
}
// from screen to world coordinate ; linear mapping
// uses global cons
double GiveWx ( int ix)
{
return (WxMin + ix * wPixelWidth);
}
// uses globaal cons
double GiveWy (int iy) {
return (WyMax - iy * wPixelHeight);
} // reverse y axis
complex double GiveW( int ix, int iy){
double Wx = GiveWx(ix);
double Wy = GiveWy(iy);
return Wx + Wy*I;
}
int SetZPlane(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;
}
// ****************** DYNAMICS = trap tests ( target sets) ****************************
// compute radius of circle around finite attractor which is independent of the image size ( iWidth/2000.0 )
// input k is a number of pixels ( in case of iWidth = 2000 )
double GiveAR(const double k){
return k*PixelWidth*iWidth/2000.0 ;
}
/* find such AR for internal LCM/J and LSM that level curves croses critical point and it's preimages
for attracting ( also weakly attracting = parabolic) dynamics
it may fail if one iteration is bigger then smallest distance between periodic point and Julia set
*/
double GiveTunedAR(int i_Max){
complex double z= zcr0; // criical point
int i;
//int i_Max = 1000;
// critical point escapes very fast here. Higher valus gives infinity
for (i=0; i< i_Max; ++i ){
z= z*z*z +c*z; // forward iteration
}
double r = cabs(z-za0);
if ( r > AR_max ) {r = AR_max;}
return r;
}
double GiveMaxFinalRadius(){
complex double z = ZxMax + ZyMax*I;
double r = log(cabs(z))/loger - 1.0; // final_z_abs in not in [0,1]
return r;
}
double GiveNormalizedFinalRadius(complex double z){
double FinalRadius = log(cabs(z))/loger - 1.0; // final_z_abs in not in [0,1]
return (FinalRadius/ MaxFinalRadius);
}
// bailout test
// z escapes when
// abs(z)> ER or cabs2(z)> ER2
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#Boolean_Escape_time
// this function is not used !!!! dead code
int Escapes(complex double z){
// here target set (trap) is the exterior circle with radsius = ER ( EscapeRadius)
// with ceter = origin z= 0
// on the Riemann sphere it is a circle with point at infinity as a center
if (cabs(z)>ER) return 1;
return 0;
}
// =====================
int IsPointInsideTrap0(complex double z){
if (cabs(z - za0) <AR) {return 1;} // circle around periodic point
return 0; // outside
}
// =====================
int IsPointInsideTrap1(complex double z){
if (cabs(z - za1) <AR) {return 1;} // circle around periodic point
return 0; // outside
}
/* ----------- 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;
}
// ***********************************************************************************************
// ********************** 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)); //
// 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 */
fprintf(stderr, "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;
}
// ============================= tests ============================================
// 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;
}
int IsInsideWWindow(complex double w){
if ( creal(w) < WxMax && creal(w) > WxMin &&
cimag(w) < WyMax && cimag(w) > WyMin) {return 1;}
return 0;
}
/*
Array A should have image of z-plane ( not w-plane)
compare of image of array A unchanged
image of w window shows part of z window and outside of z-window
"Note that the flower-shaped hole in the center is originally the edge boundary of the grid."
http://xahlee.info/SpecialPlaneCurves_dir/Inversion_dir/inversion.html
https://mathworld.wolfram.com/ConformalMapping.html
http://home.iitk.ac.in/~saiwal/engineering/complex-mappings/
*/
int ShowWWindowOnZWindow(unsigned char A[] )
{
complex double z;
//double Zx, Zy; // Z= Zx+ZY*i;
complex double w;
unsigned i; /* index of 1D array */
unsigned int ix, iy; // pixel coordinate
fprintf(stderr, "compute image ShowWWindowOnZWindow\n");
// for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy, i, w, z) 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){
z = GiveZ(ix,iy); // from screen to world coordinate
w = 1/z; // invert complex plane z
if (IsInsideWWindow(w)){
i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
A[i]=255-A[i]; // marking w window on z window
}
}
}
return 0;
}
// ------------------------------------------------------------------------------
int IsInsideZWindow(complex double z){
if ( creal(z) < ZxMax && creal(z) > ZxMin &&
cimag(z) < ZyMax && cimag(z) > ZyMin) {return 1;}
return 0;
}
/*
Array A should have image of w-plane ( not z-plane)
compare of image of array A unchanged
image of w window shows part of z window and outside of z-window
"Note that the flower-shaped hole in the center is originally the edge boundary of the grid."
http://xahlee.info/SpecialPlaneCurves_dir/Inversion_dir/inversion.html
https://mathworld.wolfram.com/ConformalMapping.html
http://home.iitk.ac.in/~saiwal/engineering/complex-mappings/
*/
int ShowZWindowOnWWindow(unsigned char A[] )
{
complex double z;
//double Zx, Zy; // Z= Zx+ZY*i;
complex double w;
unsigned i; /* index of 1D array */
unsigned int ix, iy; // pixel coordinate
fprintf(stderr, "compute image ShowZWindowOnWWindow\n");
// for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy, i, w, z) 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){
w = GiveW(ix,iy); // from screen to world coordinate
z = 1/w; // invert complex plane z
if (IsInsideZWindow(z)){
i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
A[i]=255-A[i]; // marking w window on z window
}
}
}
return 0;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// ++++++++++++++++++++++++++++++++++++++++++ color +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
unsigned char ComputeColorOfFatou (complex double z )
{
int nMax = iterMax;
int n;
double r;
for (n=0; n < nMax; n++){ //forward iteration
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
r =cabs(z);
if (r > ER) // esaping = exterior
{
//uExterior += 1;
return iColorOfExterior;
}
if ( IsPointInsideTrap0(z)) {
//uInterior +=1;
return iColorOfInterior1;}
if (IsPointInsideTrap1(z)){
//uInterior +=1;
return iColorOfInterior2;}
}
return iColorOfUnknown;
}
// ***************************************************************************************************************************
// ************************** internal LSM/J*****************************************
// ****************************************************************************************************************************
unsigned char ComputeColorOfIntLSM (complex double z)
{
double r;
int nMax = iterMax;
int n;
for (n=0; n < nMax; n++){ //forward iteration
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
r =cabs(z);
if (r > ER) // esaping = exterior
{
//uExterior += 1;
return iColorOfExterior;
}
// internal level sets the same for both sets
if ( IsPointInsideTrap0(z) || IsPointInsideTrap1(z)) {
if ( n % 2 )
{return (n*10) % 255 ;}
else {return (n*11) % 255;}
}
}
//uUnknown += 1;
return iColorOfUnknown;
}
// ***************************************************************************************************************************
// ************************** external LSM/J*****************************************
// ****************************************************************************************************************************
int GiveEscapeTime(complex double z){
int nMax = iterMax_LSM;
double cabsz;
int n;
for (n=0; n < nMax; n++){ //forward iteration
cabsz = cabs(z);
if (cabsz > ER_LSM) {break;} // escaping = exterior
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
}
return n;
}
unsigned char ComputeColorOfExtLSM(complex double z){
unsigned char iColor;
int n; // escape time
n = GiveEscapeTime(z);
// manually udjusted series of ordered colors ( shades of gray )
if (n==iterMax_LSM)
{ iColor = iColorOfInterior;}
else iColor = 255 - 230.0*((double) n)/18.0; // nMax or lower values in denominator
return iColor;
}
// ***************************************************************************************************************************
// ************************** LSM/J = both external and internal *****************************************
// ****************************************************************************************************************************
int GiveEscapeAndAttractionTime(complex double z){
int nMax = iterMax_LSM;
double cabsz;
int n;
for (n=0; n < nMax; n++){ //forward iteration
cabsz = cabs(z);
if (cabsz > ER_LSM) break; // escaping = exterior
if ( IsPointInsideTrap0(z) || IsPointInsideTrap1(z)) {break;} // attracted to finite attractor = interior
//if (cabsz< PixelWidth) break; // fails into finite attractor = interior, but not for disconnected Julia sets, then critical point and its preimages !!!!
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
}
return n;
}
unsigned char ComputeColorOfLSM(complex double z){
unsigned char iColor;
int n; // escape time
n = GiveEscapeAndAttractionTime(z);
// manually udjusted series of ordered colors ( shades of gray )
iColor = 255 - 230.0*((double) n)/18.0; // nMax or lower values in denominator
return iColor;
}
// ***************************************************************************************************************************
// ************************** DEM/J*****************************************
// ****************************************************************************************************************************
unsigned char ComputeColorOfDEM(const complex double z0){
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#DEM.2FJ
int nMax = iterMax_DEM;
complex double z = z0;
complex double dz = 1.0; // is first derivative with respect to z.
//double distance;
ER_DEM = 4.0;
double cabsz;
int n;
for (n=0; n < nMax; n++){ //forward iteration
cabsz = cabs(z);
if ( IsPointInsideTrap0(z) || IsPointInsideTrap1(z)) {return iColorOfInterior;} // attracted to finite attractor = interior
if (cabsz > ER_DEM || cabs(dz)> 1e60) {break;} // big values
//dz = (3.0*z*z + c)*dz;
z = z*z*z + c*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;
}
// ***************************************************************************************************************************
// ************************** only boundary by DEM/J*****************************************
// ****************************************************************************************************************************
unsigned char ComputeColorOfDEMJ_boundary(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 (cabsz< PixelWidth) return iColorOfInterior; // falls into finite attractor = interior
dz = (3.0*z*z + c)*dz;
z = z*z*z +c*z ; /* forward iteration : complex qubic polynomial */
}
distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
if (distance <distanceMax) return iColorOfBoundary; // distanceMax = BoundaryWidth*PixelWidth;
// else
return iColorOfExterior;
}
// plots raster point (ix,iy)
int DrawPointOfDEMJ_boundary (unsigned char A[], int PlaneInversion, int ix, int iy, unsigned char iColor0)
{
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 */
if (PlaneInversion)
{
complex double w;
w = GiveW(ix,iy);
z = 1/w;
}
else { z = GiveZ(ix,iy);}
iColor = ComputeColorOfDEMJ_boundary(z);
if (iColor == iColorOfBoundary) // check if it is boundary
{ A[i] = iColor0 ;} // draw only boundary without changing other parts using color iColor0
return 0;
}
// fill array
// uses global var : ...
// scanning complex plane
int DrawImageOfDEMJ_boundary (unsigned char A[], int PlaneInversion, const unsigned char iColor)
{
unsigned int ix, iy; // pixel coordinate
fprintf(stderr, "compute image DEM boundary\n");
// for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy) 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)
DrawPointOfDEMJ_boundary(A, PlaneInversion, ix, iy, iColor); //
}
return 0;
}
// ***************************************************************************************************************************
// ************************** Unknown: boundary and slow dynamics *****************************************
// ****************************************************************************************************************************
unsigned char ComputeColorOfUnknown(complex double z){
int nMax = 20; // very low value
double cabsz;
int n;
for (n=0; n < nMax; n++){ //forward iteration
cabsz = cabs(z);
if (cabsz > 10000000000*ER ) return iColorOfExterior; // big values
if (cabsz < (PixelWidth/100)) return iColorOfInterior; // falls into finite attractor = interior
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
}
//printf("found \n");
return iColorOfUnknown;
}
// ***************************************************************************************************************************
// ************************** binary decomposition BD/J*****************************************
// ****************************************************************************************************************************
unsigned char ComputeColorOfBD(complex double z){
int nMax = iterMax_LSM;
double cabsz;
unsigned char iColor;
int n;
for (n=0; n < nMax; n++){ //forward iteration
cabsz = cabs(z);
if (cabsz > ER_LSM) break; // esacping
//if (cabsz< PixelWidth) break; // fails into finite attractor = interior but not for disconnected Julia sets, then critical point and its preimages !!!!
// attracted to finite attractor = interior
if ( IsPointInsideTrap0(z)) { if (cimag(z)>za0im) {return 255;} else return 0;}
if ( IsPointInsideTrap1(z)) { if (cimag(z)>za1im) {return 255;} else return 0;}
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
}
// exterior
if (cimag(z)>0.0)
iColor = 255;
else iColor = 0;
return iColor;
}
// ***************************************************************************************************************************
// ************************** modified binary decomposition MBD *****************************************
// ****************************************************************************************************************************
unsigned char ComputeColorOfMBD(complex double z){
// const number of iterations
int nMax = 7;
//double cabsz;
unsigned char iColor;
int n;
for (n=0; n < nMax; n++){ //forward iteration
//cabsz = cabs(z);
//if (cabsz > ER) break; // esacping
//if (cabsz< PixelWidth) break; // falls into finite attractor = interior
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
}
//if (cabs(z) > 2.0)
{ // exterior
if (creal(z)>0.0)
iColor = 255;
else iColor = 0;
}
// else iColor = iColorOfInterior;
return iColor;
}
// ***************************************************************************************************************************
// ************************** binary decomposition boundaries with texture mapping *****************************************
// ****************************************************************************************************************************
// https://fractalforums.org/programming/11/how-many-different-ways-are-there-to-show-such-set/3874
/*
to add
https://www.iquilezles.org/www/articles/distfunctions2d/distfunctions2d.htm
2D Gray gradient = 2D gray texture
input x and y is in [0,1] range
*/
double Give2DGrayGradient(double x, double y, const int k ){
double d; // position of the color in the gradient . It is in [0,1]] range
switch(k){
case 0: {d = max(fabs(x - 0.5) ,fabs(y-0.5)); break;} //
case 1: {d = min(x,y); break;}
case 2 : {d = fabs(x)+fabs(y) -0.5; break;}
case 3 : {d = y; break;}
case 4 : {d = x; break;}
// gradients 5,6,7 are similar , difference : 1, 1,5, 2.0
case 5: {x =x - 0.5; y =y - 0.5; d = cabs(x+y*I); break;} // cabs(z)
case 6: {x =1.5*(x - 0.5); y =1.5*(y - 0.5); d = cabs(x+y*I); break;} // cabs(z)
case 7: {x =2.0*(x - 0.5); y =2.0*(y - 0.5); d = cabs(x+y*I); break;} // cabs(z)
default:{ d= 0.0; }
}
return d;
}
unsigned char ComputeColorOfTexture(complex double z, const int k){
int nMax = iterMax_LSM;
double cabsz;
unsigned char iColor;
int n;
for (n=0; n < nMax; n++){ //forward iteration
cabsz = cabs(z);
if (cabsz > ER_LSM) break; // esacping
//if (cabsz< PixelWidth) break; // fails into finite attractor = interior but not for disconnected Julia sets, then critical point and its preimages !!!!
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
}
// https://gitlab.com/adammajewski/mandelbrot-book_book/-/blob/master/README.md#final-angle
//if (n < nMax) // exterior
// {
//double et = ((double)n)/nMax; // ok but the same for all points inside level set so segmentation
double final_angle = c_turn(z); // in [0,1] range
//double final_radius = GiveNormalizedFinalRadius(z); // = final_absz should be in [0,1]
//}
double y = frac(n-log(log(cabsz)));
// inside each cell point has additional coordinate w = (final_angle, final_radius) in [0,1]x[0,1]
double gray = Give2DGrayGradient(final_angle, y, k);
iColor = gray*255;
// bd : mark each cell
//if (cimag(z)>0.0) iColor =255 -iColor;
return iColor;
}
// plots raster point (ix,iy)
int DrawPointOfTexture (unsigned char A[], int ix, int iy, const int k)
{
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 = ComputeColorOfTexture(z, k);
A[i] = iColor ; // interior
return 0;
}
// fill array
// uses global var : ...
// scanning complex plane
int DrawImageOfTexture (unsigned char A[], const int k)
{
unsigned int ix, iy; // pixel coordinate
fprintf(stderr, "compute image texture k = %d\n", k);
// for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy) 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)
DrawPointOfTexture(A, ix, iy, k); //
}
return 0;
}
// ***********************************************************************************************
//*************************************** SAC/J **************************************************
// *****************************************************************************************
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/stripeAC
// SAC = Stripe Average Coloring
//
// the addend function
// input : complex number z
// output : double number t
double Give_t(double complex z){
return 0.5+0.5*sin(s*carg(z));
}
/*
input :
- complex number
- intege
output = average
*/
double Give_Arg(double complex z , int iMax)
{
int i=0; // iteration
//double complex Z= 0.0; // initial value for iteration Z0
double A = 0.0; // A(n)
double prevA = 0.0; // A(n-1)
double R; // =radius = cabs(Z)
double d; // smooth iteration count
double complex dz = 1.0; // first derivative with respect to z
double de; // Distance Estimation from DEM/J
// iteration = computing the orbit
for(i=0;i<iMax;i++)
{
dz = (3.0*z*z + c)*dz;
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
if (i>i_skip) A += Give_t(z); //
R = cabs(z);
// if(R > EscapeRadius) break; // exterior of M set
if (R > 1e60 || cabs(dz)> 1e60) break; // prevent NAN
prevA = A; // save value for interpolation
} // for(i=0
if (i == iMax)
A = -1.0; // interior
else { // exterior
de = 2 * R * log(R) / cabs(dz);
if (de < distanceMax) A = FP_ZERO; // boundary
else {
// computing interpolated average
A /= (i - i_skip) ; // A(n)
prevA /= (i - i_skip - 1) ; // A(n-1)
// smooth iteration count
d = i + 1 + log(lnER/log(R))/M_LN2;
d = d - (int)d; // only fractional part = interpolation coefficient
// linear interpolation
A = d*A + (1.0-d)*prevA;
}
}
return A;
}
unsigned char ComputeColorOfSAC(complex double z){
unsigned char iColor;
double arg;
arg = Give_Arg( z, 2500); // N in wiki
// color is proportional to arg
if (arg < 0.0)
iColor = 0; // interior
else //
{if (arg == FP_ZERO)
iColor = 255; // boundary
else iColor = (unsigned char) (255 - 255*arg );// exterior
}
return iColor;
}
// ***************************************************************************************************************************
// ************************** DLD/J*****************************************
// ****************************************************************************************************************************
/* partial pnorm
input: z , zn = f(z), p
output ppn
*/
double
ppnorm (complex double z, complex double zn, double p)
{
double s[2][3]; // array for 2 points on the Riemann sphere
int j;
double d; // denominator
double x;
double y;
double ds;
double ppn = 0.0;
// map from complex plane to riemann sphere
// z
x = creal (z);
y = cimag (z);
d = x * x + y * y + 1.0;
s[0][0] = (2.0 * x) / d;
s[0][1] = (2.0 * y) / d;
s[0][2] = (d - 2.0) / d; // (x^2 + y^2 - 1)/d
// zn
x = creal (zn);
y = cimag (zn);
d = x * x + y * y + 1.0;
s[1][0] = (2.0 * x) / d;
s[1][1] = (2.0 * y) / d;
s[1][2] = (d - 2.0) / d; // (x^2 + y^2 - 1)/d
// sum
for (j = 0; j < 3; ++j)
{
ds = fabs (s[1][j] - s[0][j]);
// normal: neither zero, subnormal, infinite, nor NaN
//if (fpclassify (ds) !=FP_INFINITE)
//if (isnormal(ds))
// it is solved by if (cabs(z) > 1e60 ) break; procedure in parent function
ppn += pow (ds, p); // |ds|^p
// else {ppn = 10000.0; printf("ds = infty\t");} //
}
return ppn;
}
// DLD = Discret Lagrangian Descriptior
double
lagrangian (complex double z0, complex double c, int iMax, double p)
{
int i; // number of iteration
double d = 0.0; // DLD = sum
double ppn; // partial pnorm
complex double z = z0;
complex double zn; // next z
for (i = 0; i < iMax; ++i)
{
zn = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
ppn = ppnorm (z, zn, p);
d += ppn; // sum
//
z = zn;
//if (! isnormal(d)) { return 0.0; } // not works
if (cabs (z) > ER_DLD ) //1e6)
break; // exterior : big values produces artifacts on the image
}
//if (d<0.0) {// interior
// d(z1a) - d(z21) = -0.0804163521959989
// d = - d;
// d = (db - d) /dd ; // normalize, see test_interior
//d = d*d;
//if (d>1.0) {printf("d int > 1.0\n");
/// }
// else {
d = d / ((double) i); // averaging not summation
//d = d*me;} // exterior
return d;
}
unsigned char
ComputeColorOfDLD (complex double z)
{
//double cabsz;
int iColor;
double d;
int N = iterMax_DLD; // N in wiki = fixed number : maximal number of iterations
//if (FatouType == 1)
// { // interior
// d = lagrangian (z, c, N, p);
// modify gradient position
//{d = d - (int)d;} // only fractional part
// d = d * d * mi;
//if ( d< 1.0 ) d = 0.0;
// } //
//else
//{
d = lagrangian (z, c, N, p); //
//}
iColor = (int) (d * 255) % 255; // nMax or lower walues in denominator
return (unsigned char) iColor;
}
//=========================================
// ***************************************************************************************************************************
// ************************** NPM/J = Normal Potential *****************************************
// ****************************************************************************************************************************
/*
The dot product of two vectors a = [a1, a2, ..., an] and b = [b1, b2, ..., bn] is defined as:[1]
d = a1b1 + a2b2
*/
double cdot(double complex a, double complex b){
return creal(a)*creal(b)+cimag(a)*cimag(b);
}
//
// output
//
double GiveReflection(double complex z )
{
int i=0; // iteration
int iMax = 2000;
// https://en.wikibooks.org/wiki/Fractals/Mathematics/Derivative
double complex dz = 1.0; // derivative with respect to z
double reflection = 0.0; //
double h2 = 1.5 ; // height factor of the incoming light
double angle = 45.0/360.0 ; // incoming direction of light in turns
double complex v = cexp(2.0*angle *M_PI* I); // = exp(1j*angle*2*pi/360) // unit 2D vector in this direction
// incoming light 3D vector = (v.re,v.im,h2)
// https://en.wikipedia.org/wiki/Lambertian_reflectance
double complex u;
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
dz = 1.0;
for(i=0;i<iMax;i++)
{
dz = (3.0*z*z + c)*dz;
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
if(cabs(z) > ER_NP)
{ // exterior
u = z / dz;
u = u / cabs(u);
reflection = cdot(u, v) + h2; /* use the simplest model for the shading:
Lambert, which consists in using the dot product of (x,y,1) with a constant vector indicating the direction of the light. */
reflection = reflection/(1.0 + h2); // rescale so that t does not get bigger than 1
if (reflection<0.0) reflection =0.0;
break;
}
}
return reflection;
}
// Potential to color
unsigned char ComputeColorOfNP(complex double z){
//https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set#Normal_map_effect
double reflection;
unsigned char iColor;
// compute
reflection = GiveReflection( z);
//
//if (reflection < )
//{ /* interior */
// iColor = 0;}
//else // exterior
{ iColor = 255 * reflection;}
return iColor;
}
// https://en.wikipedia.org/wiki/Shading
// normal = perpendicular
// shading using Normal map and Potential
// https://en.wikipedia.org/wiki/Lambertian_reflectance
// http://www.math.titech.ac.jp/~kawahira/gallery/movies/movies.html
// see 0_1.avi and image
//
// ***************************************************************************************************************************
// ************************** NDM/J = Normal Distance *****************************************
// ****************************************************************************************************************************
// normal = perpendicular
// shading using Normal map and Potential
// https://en.wikipedia.org/wiki/Lambertian_reflectance
// http://www.math.titech.ac.jp/~kawahira/gallery/movies/movies.html
// see 0_1.avi and image
//
//
// output
//
double GiveReflectionD(double complex z )
{
int i=0; // iteration
int iMax = 2000;
// https://en.wikibooks.org/wiki/Fractals/Mathematics/Derivative
double complex dz = 1.0; // first derivative with respect to z
double complex dz2 = 0.0; // second derivative with respect to z
double reflection = 0.0; //
double lo;
double h2 = 1.5 ; // height factor of the incoming light
double angle = 45.0/360.0 ; // incoming direction of light in turns
double complex v = cexp(2.0*angle *M_PI* I); // = exp(1j*angle*2*pi/360) // unit 2D vector in this direction
// incoming light 3D vector = (v.re,v.im,h2)
// https://en.wikipedia.org/wiki/Lambertian_reflectance
double complex u;
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
dz = 1.0;
dz2 = 0.0;
for(i=0;i<iMax;i++)
{
dz2 = 2.0* ( dz2*z + dz*dz);//2*(der2*z+der**2)
dz = (3.0*z*z + c)*dz;
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
if(cabs(z) > ER_NP)
{ // exterior
/*
lo = 0.5*log(squared_modulus(z))
u = z*der*((1+lo)*conj(der**2)-lo*conj(z*der2))
u = u/abs(u)
*/
lo = 0.5*log(cabs(z));
u = z*dz*((1.0+lo)*conj(dz*dz)-lo*conj(z*dz2));
//u = z / dz;
u = u / cabs(u);
reflection = cdot(u, v) + h2; // use the simplest model for the shading: Lambert, which consists in using the dot product of (x,y,1) with a constant vector indicating the direction of the light.
reflection = reflection/(1.0 + h2); // rescale so that t does not get bigger than 1
if (reflection<0.0) reflection =0.0;
break;
}
}
return reflection;
}
// Distance to color
unsigned char ComputeColorOfND(complex double z){
//https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set#Variation
double reflection;
unsigned char iColor;
// compute
reflection = GiveReflectionD( z);
//
//if (reflection < )
//{ /* interior */
// iColor = 0;}
//else // exterior
{ iColor = 255 * reflection;}
return iColor;
}
// -------------------------- potential========
double ComputePotential(const complex double z0){
double potential = 0.0; // interior
double s = 0.5;
complex double z = z0;
double r;
int iter;
for (iter = 0; iter < iterMax_pot; ++iter){
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
s *= 0.5; //
r = cabs(z);
if (r > ER_pot) {break;}
}
potential = s*log2(r); // log(zn)* 2^(-n)
return potential;
}
unsigned char ComputeColorOfPOT(complex double z){
double potential = ComputePotential(z);
if (PlaneInversion) // usung global var
{potential /= 4.0;} // manual normalize
unsigned char iColor = 255 * sqrt(sqrt(potential));
return iColor;
}
double GiveSmoothEscapeTime(complex double z){
int nMax = iterMax_LSM;
double cabsz;
int n;
for (n=0; n < nMax; n++){ //forward iteration
cabsz = cabs(z);
if (cabsz > ER_LSM) break; // esacping
//if (cabsz< PixelWidth) break; // fails into finite attractor = interior, but not for disconnected Julia sets, then critical point and its preimages !!!!
z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */
}
//float sn = n - log2(log2(dot(z,z))) + 4.0; // equivalent optimized smooth iteration count
double sn = ER_LSM/cabsz;
//n- log2(log2(cdot(z,z))) + 4.0;
//sn = sn / nMax; // map to [0,1] range
return sn;
}
//
unsigned char ComputeColorOfBlend(complex double z){
double SET = GiveSmoothEscapeTime(z);
SET = sqrt(SET);
SET = 1.0 - SET;
//
double ColorSET = SET*255;
//
double ColorNP = ComputeColorOfNP(z);
unsigned char iColor = (ColorSET+ ColorNP)/ 2.0; // average blend mode
return iColor;
}
/*
int local_setup(int PlaneInversion){
if (PlaneInversion)
{ MaxImagePotential =ComputePotential( 1.0/ 0.0);}
//else {MaxImagePotential}
return 0;
}
*/
/* ==================================================================================================
============================= 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 = ComputeColorOfDEM(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 (FunctionTypeT FunctionType, int PlaneInversion, unsigned char A[], 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 */
if (PlaneInversion)
{
complex double w;
w = GiveW(ix,iy);
z = 1/w;
}
else { z = GiveZ(ix,iy);}
iColor = ComputeColor(FunctionType, z);
A[i] = iColor ; //
return 0;
}
// fill array
// uses global var : ...
// scanning complex plane
int DrawImage (FunctionTypeT FunctionType, int PlaneInversion, unsigned char A[])
{
unsigned int ix, iy; // pixel coordinate
//local_setup(PlaneInversion);
fprintf(stderr, "compute image FunctionType = %d PlaneInversion = %d\n", FunctionType, PlaneInversion);
// for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy) 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)
DrawPoint(FunctionType, PlaneInversion, A, ix, iy); //
}
return 0;
}
// ------------------------------------------------------------------------------------------------------------
// ---------------------------Additional drawings : traps, attractors and critical orbits -------------------------
// --------------------------------------------------------------------------------------------------------------
int PointIsInsideZRectangle(const double complex z){
if (ZxMin < creal(z) && creal(z)<ZxMax && ZyMin < cimag(z) && cimag(z)<ZyMax)
{return 1;}
return 0;
}
int IsInside (int x, int y, int xcenter, int ycenter, int r){
double dx = x- xcenter;
double dy = y - ycenter;
double d = sqrt(dx*dx+dy*dy);
if (d<r)
return 1;
return 0;
}
int PlotBigPoint(complex double z, unsigned char A[]){
unsigned int ix_seed = (creal(z)-ZxMin)/PixelWidth;
unsigned int iy_seed = (ZyMax - cimag(z))/PixelHeight;
unsigned int i;
/* mark seed point by big pixel */
int iSide =3.0*iWidth/2000.0 ; /* half of width or height of big pixel */
int iY;
int iX;
for(iY=iy_seed-iSide;iY<=iy_seed+iSide;++iY){
for(iX=ix_seed-iSide;iX<=ix_seed+iSide;++iX){
if (IsInside(iX, iY, ix_seed, iy_seed, iSide)) {
i= Give_i(iX,iY); /* index of _data array */
//A[i]= 255-A[i]; // inverted color
A[i] = 0; // allways black
}}}
return 0;
}
// =====================
int IsPointInsideTraps(unsigned int ix, unsigned int iy){
complex double z = GiveZ (ix, iy);
if ( IsPointInsideTrap0(z)) {return 1;} // circle with prabolic point on it's boundary
if (IsPointInsideTrap1(z)) {return 1;}
return 0; // outside
}
int MarkTraps(unsigned char A[]){
unsigned int ix, iy; // pixel coordinate
unsigned int i;
fprintf (stderr, "Mark traps \n");
// for all pixels of image
#pragma omp parallel for schedule(dynamic) private(ix,iy) 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){
if (IsPointInsideTraps(ix, iy)) {
i= Give_i(ix,iy); /* index of _data array */
A[i]= 255-A[i]; // inverse color
}}}
return 0;
}
int DrawForwardOrbit(const complex double z0, const unsigned long long int iMax, unsigned char A[] )
{
unsigned long long int i; /* nr of point of critical orbit */
complex double z = z0;
PlotBigPoint(z, A);
/* forward orbit of critical point */
for (i=1;i<iMax ; ++i)
{
z = z*z*z+ c*z;
// if (cabs2(z - z2a) > 2.0) {return 1;} // escaping
if (PointIsInsideZRectangle(z))
{PlotBigPoint(z, A);}
else fprintf(stderr, "bad point z\n");
}
return 0;
}
void DrawCriticalOrbits(unsigned char A[]){
unsigned long long int iMax = iterMax_LSM ;
DrawForwardOrbit(zcr0, iMax, A);
DrawForwardOrbit(zcr1, iMax, A);
}
//-----------------------------------------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------------------------------------
// ------------------------------------------------------------------------------------------------------------------------------
int Test()
{
unsigned int ix, iy; // pixel coordinate
complex double z;
double SET;
int ET;
//local_setup(PlaneInversion);
fprintf(stderr, "test\n");
// for all pixels of image
ix = 0;
for (iy = iyMin; iy <= iyMax; ++iy){
z = GiveZ(ix,iy);
ET = GiveEscapeTime(z);
SET = GiveSmoothEscapeTime(z);
printf(" %d \t %d \t %f \n", iy, ET, SET);
ix = ix +1;
//
}
return 0;
}
// *******************************************************************************************
// ********************************** save A array to pgm file ****************************
// *********************************************************************************************
int
SaveArray2PGMFile (unsigned char A[], char *shortName , char *comment)
{
FILE *fp;
const unsigned int MaxColorComponentValue = 255; /* color component is coded from 0 to 255 ; it is 8 bit color file */
// https://programmerfish.com/create-output-file-names-using-a-variable-in-c-c/
char fileName[512];
const char* fileType = ".pgm";
sprintf(fileName,"%s%s", shortName, fileType); //
char long_comment[200];
sprintf (long_comment, "f(z) = z*z*z +c*z where c = %f %+f*i ; %s", creal(c), cimag(c),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 %u requested\n", rSize, iSize);}
//
NumberOfImages +=1; // count images using global variable
return 0;
}
int PrintInfoAboutProgam()
{
printf("Number of pgm images = %d \n", NumberOfImages);
// display info messages
printf ("Numerical approximation of Julia set for fc(z)= z^2 + c \n");
//printf ("iPeriodParent = %d \n", iPeriodParent);
//printf ("iPeriodOfChild = %d \n", iPeriodChild);
printf ("parameter c = %.16f %+.16f*i \n", creal(c), cimag(c));
printf ("Image Width = %f in world coordinate\n", ZxMax - ZxMin);
printf ("PixelWidth = %f \n", PixelWidth);
printf("for DEM/J \n");
if ( distanceMax<0.0 || distanceMax > ER ) printf("bad distanceMax\n");
printf("Max distance from exterior to the boundary = distanceMax = %.16f = %f pixels\n", distanceMax, BoundaryWidth);
printf("\n");
// image corners in world coordinate
// center and radius
// center and zoom
// GradientRepetition
printf ("Maximal number of iterations = iterMax = %ld \n", iterMax);
printf ("For LSM/J \n");
printf ("Maximal number of iterations = iterMax_LSM = %ld \n", iterMax_LSM);
printf ("Escape Radius = ER_LSM = %f \n", ER_LSM);
printf("\n");
printf ("ratio of image = %f ; it should be 1.000 ...\n", ratio);
//
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
return 0;
}
// *****************************************************************************
//;;;;;;;;;;;;;;;;;;;;;; setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
// **************************************************************************************
int setup ()
{
fprintf (stderr, "setup start\n");
// c= 1.0049542069308062 +0.1008317508132964*i
c = 1.02*cexp(0.1*I); ; // https://web.archive.org/web/20161024194536/http://www.ijon.de/mathe/julia/some_julia_sets_3.html
/* 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].
center = za1;
SetZPlane( 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 ...
wPixelWidth = (WxMax-WxMin)/ixMax;
wPixelHeight =(WyMax-WyMin)/iyMax;
//ER2 = ER * ER; // for numerical optimisation in iteration
lnER = log(EscapeRadius); // ln(ER)
loger = log(ER_LSM); // for texture
ER_LSM = 3.0; //GiveER(10); // find such ER for LSM/J that level curves croses critical point and it's preimages
ER_DLD = 3.0; //GiveER(7);
MaxFinalRadius = GiveMaxFinalRadius();
AR_max = GiveAR(50.0); //PixelWidth*50.0*iWidth/2000.0 ; // moved to main
// BD
za0im = cimag(za0);
za1im = cimag(za1);
/* create dynamic 1D arrays for colors ( shades of gray ) */
data = malloc (iSize * sizeof (unsigned char));
edge = malloc (iSize * sizeof (unsigned char));
edge2 = malloc (iSize * sizeof (unsigned char));
//
if (data == NULL || edge == NULL || edge2 == NULL ){
fprintf (stderr, " Could not allocate memory");
return 1;
}
BoundaryWidth = 1.0*iWidth/2000.0 ; // measured in pixels ( when iWidth = 2000)
distanceMax = BoundaryWidth*PixelWidth;
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);
free(edge2);
PrintInfoAboutProgam();
return 0;
}
// ********************************************************************************************************************
/* ----------------------------------------- main -------------------------------------------------------------*/
// ********************************************************************************************************************
int main () {
setup ();
PlaneInversion = 0;
/*
printf("bigger trap radius for speed \n");
AR = AR_max;
printf("AR_max = %.16f = %f pixels\n",AR_max , AR_max/PixelWidth);
DrawImage(Fatou, PlaneInversion, data);
SaveArray2PGMFile (data, "Fatou", "Fatou set");
ComputeBoundaries(data,edge);
SaveArray2PGMFile (edge, "FatouB", "boundaries of the Fatou set ");
CopyBoundaries(edge,data);
SaveArray2PGMFile (data, "FatouAndB", "Fatou set and boundaries");
// mark attractors
PlotBigPoint(za0, data);
PlotBigPoint(za1, data);
SaveArray2PGMFile (data, "Fatou_p", "Fatou and attracting points");
MarkTraps(data);
SaveArray2PGMFile (data, "Fatou_traps", "Fatou and traps");
DrawCriticalOrbits(data);
SaveArray2PGMFile (data, "Fatou_cr", "Fatou, traps and critical orbits");
AR = GiveAR(10);
DrawImage(BD, PlaneInversion, data);
SaveArray2PGMFile (data, "BD", "BD"); // name of the file is name.png
*/
printf("smaller trap radius for more detailes \n");
AR = GiveTunedAR(400);
printf("tuned AR = %.16f = %.16f *AR_max = %f pixels\n",AR, AR/AR_max, AR/PixelWidth);
/*
DrawImage(IntLSM, PlaneInversion, data);
SaveArray2PGMFile (data, "IntLSM", "internal level sets = level sets of attraction time"); // name of the file is name.png
ComputeBoundaries(data,edge);
SaveArray2PGMFile (edge, "IntLCM", "internal Level Curves = boundaries of Int Level Sests "); // name of the file is name.png
CopyBoundaries(edge,data);
SaveArray2PGMFile (data, "IntLSC", "Int LevelSets and it's boundaries = LevelCurves"); // name of the file is name.png
*/
DrawImage(LSM, PlaneInversion, data);
SaveArray2PGMFile (data, "LSM", "internal and external level sets = level sets of attraction and escaping time"); // name of the file is name.png
ComputeBoundaries(data,edge);
SaveArray2PGMFile (edge, "LCM_zoom", "internal and external Level Curves = boundaries of both Level Sests "); // name of the file is name.png
DrawCriticalOrbits(edge);
SaveArray2PGMFile (edge, "LCM_cr_zoom", "LCM and critical orbits");
/*
CopyBoundaries(edge,data);
SaveArray2PGMFile (data, "LSCM", "Int and ext LevelSests and it's boundaries = LevelCurves"); // name of the file is name.png
DrawImage(DEM, PlaneInversion, data);
SaveArray2PGMFile (data, "DEM", "external Distance Estimation Method = boundary of filled Julia set"); // name of the file is name.png
*/
//Test();
//
end();
return 0;
}
Makefile
[edit]%3D_z%5E3%2B(1.0149042485835864102%2B0.10183008497976470119i)*z;_(zoom).png&action=edit§ion=4){kind=link}
all:
chmod +x d.sh
./d.sh
bash source code
[edit]%3D_z%5E3%2B(1.0149042485835864102%2B0.10183008497976470119i)*z;_(zoom).png&action=edit§ion=5){kind=link}
#!/bin/bash
# script file for BASH
# which bash
# save this file as d.sh
# chmod +x d.sh
# ./d.sh
# checked in https://www.shellcheck.net/
# display OMP info
export OMP_DISPLAY_ENV="TRUE"
printf "make pgm and txt files \n"
gcc d.c -lm -Wall -march=native -fopenmp
if [ $? -ne 0 ]
then
echo ERROR: compilation failed !!!!!!
exit 1
fi
time ./a.out > a.txt
export OMP_DISPLAY_ENV="FALSE"
printf "convert all pgm files to png using Image Magic convert \n"
# for all pgm files in this directory
for file in *.pgm ; do
# b is name of file without extension
b=$(basename "$file" .pgm)
# convert using ImageMagic
convert "${b}".pgm -resize 2000x2000 "${b}".png
echo "$file"
done
export OMP_DISPLAY_ENV="TRUE"
printf "display OMP info \n"
printf "delete all pgm files \n"
rm ./*.pgm
echo OK
# end
text output
[edit]%3D_z%5E3%2B(1.0149042485835864102%2B0.10183008497976470119i)*z;_(zoom).png&action=edit§ion=6){kind=link}
chmod +x d.sh ./d.sh make pgm and txt files OPENMP DISPLAY ENVIRONMENT BEGIN _OPENMP = '201511' OMP_DYNAMIC = 'FALSE' OMP_NESTED = 'FALSE' OMP_NUM_THREADS = '8' OMP_SCHEDULE = 'DYNAMIC' OMP_PROC_BIND = 'FALSE' OMP_PLACES = '' OMP_STACKSIZE = '0' OMP_WAIT_POLICY = 'PASSIVE' OMP_THREAD_LIMIT = '4294967295' OMP_MAX_ACTIVE_LEVELS = '2147483647' OMP_CANCELLATION = 'FALSE' OMP_DEFAULT_DEVICE = '0' OMP_MAX_TASK_PRIORITY = '0' OMP_DISPLAY_AFFINITY = 'FALSE' OMP_AFFINITY_FORMAT = 'level %L thread %i affinity %A' OPENMP DISPLAY ENVIRONMENT END setup start end of setup compute image FunctionType = 3 PlaneInversion = 0 bad point z ... bad point z allways free memory (deallocate ) to avoid memory leaks real 19m13,945s user 138m33,613s sys 0m5,662s convert all pgm files to png using Image Magic convert LCM_cr_zoom.pgm LCM_zoom.pgm LSM.pgm display OMP info delete all pgm files OK smaller trap radius for more detailes tuned AR = 0.0038665412903206 = 0.1546539185302452 *AR_max = 77.326959 pixels File LSM.pgm saved . Comment = f(z) = z*z*z +c*z where c = 1.014904 +0.101830*i ; internal and external level sets = level sets of attraction and escaping time File LCM_zoom.pgm saved . Comment = f(z) = z*z*z +c*z where c = 1.014904 +0.101830*i ; internal and external Level Curves = boundaries of both Level Sests File LCM_cr_zoom.pgm saved . Comment = f(z) = z*z*z +c*z where c = 1.014904 +0.101830*i ; LCM and critical orbits Number of pgm images = 3 Numerical approximation of Julia set for fc(z)= z^2 + c parameter c = 1.0149042485835864 +0.1018300849797647*i Image Width = 1.000000 in world coordinate PixelWidth = 0.000050 for DEM/J Max distance from exterior to the boundary = distanceMax = 0.0005000250012501 = 10.000000 pixels Maximal number of iterations = iterMax = 1000000 For LSM/J Maximal number of iterations = iterMax_LSM = 1000 Escape Radius = ER_LSM = 3.000000 ratio of image = 1.000000 ; it should be 1.000 ... gcc version: 10.2.0
references
[edit]%3D_z%5E3%2B(1.0149042485835864102%2B0.10183008497976470119i)*z;_(zoom).png&action=edit§ion=7){kind=link}
File history
Click on a date/time to view the file as it appeared at that time.
| Date/Time | Thumbnail | Dimensions | User | Comment | |
|---|---|---|---|---|---|
| current | 17:35, 17 May 2021 | | 2,000 × 2,000 (1.65 MB) | Soul windsurfer (talk | contribs) | Uploaded own work with UploadWizard |
You cannot overwrite this file.
File usage on Commons
The following 3 pages use this file:
_%3D_z*z%2B0.25.svg){kind=link}
File usage on other wikis
The following other wikis use this file:
- Usage on en.wikipedia.org
- Usage on en.wikibooks.org
%3D_z%5E3%2B(1.0149042485835864102%2B0.10183008497976470119i)*z;_(zoom).png&oldid=616389492){kind=link}