File:External rays and critical orbit landing on the parabolic fixed point for t=5 over 11.png
Page contents not supported in other languages.
Tools
Actions
General
In other projects
Appearance
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: 27 KB, MIME type: image/png)
This is a file from the Wikimedia Commons. Information from its description page there is shown below. Commons is a freely licensed media file repository. You can help. |
Summary
DescriptionExternal rays and critical orbit landing on the parabolic fixed point for t=5 over 11.png |
English: External rays and critical orbit landing on the parabolic fixed point for internal angle t=5 /11. Computed parts are white ( 255). Approximated parts are gray ( 155). Function f(z) = z*z -0.690059870015044 +0.276026482784614*I |
Date | |
Source | Own work |
Author | Soul windsurfer |
Licensing
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
/*
Adam Majewski
adammaj1 aaattt o2 dot pl // o like oxygen not 0 like zero
console program in c programing language
==============================================
---------------------------------
indent e.c
default is gnu style
-------------------
c console progam
gcc e.c -lm -Wall -march=native
time ./a.out > b.txt
gcc d.c -lm -Wall -march=native
time ./a.out
time ./a.out >e.txt
----------------------
real 0m19,809s
user 2m26,763s
sys 0m0,161s
gnuplot> set xlabel "distance"
gnuplot> set ylabel "distance"
gnuplot> set xlabel "turn"
gnuplot> set output "data.png"
gnuplot> plot "data.txt"
gnuplot> set title "distance from ray to the landing point after n iteration"
gnuplot> set output "data.png"
gnuplot> plot "data.txt"
gnuplot> unset key
gnuplot> set output "data.png"
gnuplot> plot "data.txt"
*/
#include <stdio.h>
#include <stdlib.h> // malloc
#include <string.h> // strcat
#include <stdbool.h>
#include <math.h> // M_PI; needs -lm also
#include <complex.h>
/* --------------------------------- global variables and consts ------------------------------------------------------------ */
// 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 = 2000; // proportional to size of image file !!!!!!
// The size of array has to be a positive constant integer
static unsigned 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
double ZxMin ; //= -1.7; //-2.0; //-0.05;
double ZxMax ; //= 1.7; // 2.0; //0.75;
double ZyMin ;// //= -1.7; //-2.0; //-0.1;
double ZyMax ; // = 1.7; //2.0; //0.7;
double ImageWidth;
double plane_radius = 2.0 ;
double complex plane_center = 0.0;
static double PixelWidth; // =(ZxMax-ZxMin)/ixMax;
static double PixelHeight; // =(ZyMax-ZyMin)/iyMax;
static double ratio;
// complex numbers of parametr plane
double complex c; // parameter of function fc(z)=z^2 + c
double complex a; // alfa fixed point
int child_period = 11;
int t_numerator = 5;
int t_denominator; // = child period
int e_angle_numerator0 = 341;
int e_angle_denominator = 2047;
// for external rays
int bCounted = 0; // boolean
double MinDistanceToLandingPoint = 1.0; // initial value
//int Period = 2;
unsigned long int iterMax_ray = 2000000000000; //it is also used for file name so it shlould be different values !!
unsigned long int iterMax_cr_orbit = 1000000000000; //
/* colors = shades of gray from 0 to 255 */
unsigned char iColorOfExterior = 245;
unsigned char iColorOfInterior = 55;
unsigned char iColorOfInterior1 = 210;
unsigned char iColorOfInterior2 = 180;
unsigned char iColorOfBoundary = 0;
unsigned char iColorOfUnknown = 30;
/* ------------------------------------------ functions -------------------------------------------------------------*/
int SetPlane(complex double center, double radius, double AspectRatio){
ZxMin = creal(center) - radius*AspectRatio; //-2.0; //-0.05;
ZxMax = creal(center) + radius*AspectRatio; // 2.0; //0.75;
ZyMin = cimag(center) - radius; //-2.0; //-0.1;
ZyMax = cimag(center) + radius; //2.0; //0.7;
ImageWidth = ZxMax - ZxMin;
return 0;
}
/* ----------- 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;
}
/*
gives position ( index) in 1D virtual array of 2D point Z
without bounds check !!
*/
int Give_i_from_d(complex double Z){ // double version of Give_k
/* translate from world to screen coordinate */
// iY=(ZyMax-Zy)/PixelHeight; /* */
int ix=(creal(Z)-ZxMin)/PixelWidth;
int iy=(ZyMax - cimag(Z))/PixelHeight; /* reverse Y axis */
return Give_i(ix,iy);
}
void ColorPixel(int iColor, int i, unsigned char A[])
{
A[i] = iColor;
}
int ColorPixel_d(complex double z, int iColor, unsigned char A[]){
int i = Give_i_from_d(z); // compute index of 1D array
if ( i<0 || i>iSize) {fprintf(stderr, " bad i from color pixel\n");return 1;}
ColorPixel(iColor, i, A);
//printf("plot z = %f;%f ; i = %d\n",creal(z), cimag(z), i);
return 0;
}
// plots raster point (ix,iy)
int iDrawPoint(unsigned int ix, unsigned int iy, unsigned char iColor, unsigned char A[])
{
/* i = Give_i(ix,iy) compute index of 1D array from indices of 2D array */
A[Give_i(ix,iy)] = iColor;
return 0;
}
// draws point to memmory array data
// uses complex type so #include <complex.h> and -lm
int dDrawPoint(complex double point,unsigned char iColor, unsigned char A[] )
{
unsigned int ix, iy; // screen coordinate = indices of virtual 2D array
//unsigned int i; // index of 1D array
ix = (creal(point)- ZxMin)/PixelWidth;
iy = (ZyMax - cimag(point))/PixelHeight; // inverse Y axis
iDrawPoint(ix, iy, iColor, A);
return 0;
}
/*
http://rosettacode.org/wiki/Bitmap/Bresenham%27s_line_algorithm
Instead of swaps in the initialisation use error calculation for both directions x and y simultaneously:
*/
void iDrawLine( int x0, int y0, int x1, int y1, unsigned char iColor, unsigned char A[])
{
int x=x0; int y=y0;
int dx = abs(x1-x0), sx = x0<x1 ? 1 : -1;
int dy = abs(y1-y0), sy = y0<y1 ? 1 : -1;
int err = (dx>dy ? dx : -dy)/2, e2;
for(;;){
iDrawPoint(x, y, iColor, A);
if (x==x1 && y==y1) break;
e2 = err;
if (e2 >-dx) { err -= dy; x += sx; }
if (e2 < dy) { err += dx; y += sy; }
}
}
int dDrawLineSegment(double Zx0, double Zy0, double Zx1, double Zy1, int color, unsigned char *array)
{
unsigned int ix0, iy0; // screen coordinate = indices of virtual 2D array
unsigned int ix1, iy1; // screen coordinate = indices of virtual 2D array
// first step of clipping
//if ( Zx0 < ZxMax && Zx0 > ZxMin && Zy0 > ZyMin && Zy0 <ZyMax
// && Zx1 < ZxMax && Zx1 > ZxMin && Zy1 > ZyMin && Zy1 <ZyMax )
ix0= (Zx0- ZxMin)/PixelWidth;
iy0 = (ZyMax - Zy0)/PixelHeight; // inverse Y axis
ix1= (Zx1- ZxMin)/PixelWidth;
iy1= (ZyMax - Zy1)/PixelHeight; // inverse Y axis
// second step of clipping
if (ix0 >=ixMin && ix0<=ixMax && ix0 >=ixMin && ix0<=ixMax && iy0 >=iyMin && iy0<=iyMax
&& iy1 >=iyMin && iy1<=iyMax )
iDrawLine(ix0,iy0,ix1,iy1,color, array) ;
return 0;
}
// -------------------------------
// https://en.wikipedia.org/wiki/Cohen%E2%80%93Sutherland_algorithm
typedef int OutCode;
const int INSIDE = 0; // 0000
const int LEFT = 1; // 0001
const int RIGHT = 2; // 0010
const int BOTTOM = 4; // 0100
const int TOP = 8; // 1000
// Compute the bit code for a point (x, y) using the clip rectangle
// bounded diagonally by (ZxMin, ZyMin), and (ZxMax, ZyMax)
// ASSUME THAT ZxMax, ZxMin, ZyMax and ZyMin are global constants.
OutCode ComputeOutCode(double x, double y)
{
OutCode code;
code = INSIDE; // initialised as being inside of [[clip window]]
if (x < ZxMin) // to the left of clip window
code |= LEFT;
else if (x > ZxMax) // to the right of clip window
code |= RIGHT;
if (y < ZyMin) // below the clip window
code |= BOTTOM;
else if (y > ZyMax) // above the clip window
code |= TOP;
return code;
}
// Cohen–Sutherland clipping algorithm clips a line from
// P0 = (x0, y0) to P1 = (x1, y1) against a rectangle with
// diagonal from (xmin, ymin) to (xmax, ymax).
// CohenSutherlandLineClipAndDraw
void dDrawLine(double x0, double y0, double x1, double y1,unsigned char iColor, unsigned char A[])
{
// compute outcodes for P0, P1, and whatever point lies outside the clip rectangle
OutCode outcode0 = ComputeOutCode(x0, y0);
OutCode outcode1 = ComputeOutCode(x1, y1);
bool accept = false;
while (true) {
if (!(outcode0 | outcode1)) { // Bitwise OR is 0. Trivially accept and get out of loop
accept = true;
break;
} else if (outcode0 & outcode1) { // Bitwise AND is not 0. (implies both end points are in the same region outside the window). Reject and get out of loop
break;
} else {
// failed both tests, so calculate the line segment to clip
// from an outside point to an intersection with clip edge
double x, y;
// At least one endpoint is outside the clip rectangle; pick it.
OutCode outcodeOut = outcode0 ? outcode0 : outcode1;
// Now find the intersection point;
// use formulas y = y0 + slope * (x - x0), x = x0 + (1 / slope) * (y - y0)
if (outcodeOut & TOP) { // point is above the clip rectangle
x = x0 + (x1 - x0) * (ZyMax - y0) / (y1 - y0);
y = ZyMax;
} else if (outcodeOut & BOTTOM) { // point is below the clip rectangle
x = x0 + (x1 - x0) * (ZyMin - y0) / (y1 - y0);
y = ZyMin;
} else if (outcodeOut & RIGHT) { // point is to the right of clip rectangle
y = y0 + (y1 - y0) * (ZxMax - x0) / (x1 - x0);
x = ZxMax;
} else if (outcodeOut & LEFT) { // point is to the left of clip rectangle
y = y0 + (y1 - y0) * (ZxMin - x0) / (x1 - x0);
x = ZxMin;
}
// Now we move outside point to intersection point to clip
// and get ready for next pass.
if (outcodeOut == outcode0) {
x0 = x;
y0 = y;
outcode0 = ComputeOutCode(x0, y0);
} else {
x1 = x;
y1 = y;
outcode1 = ComputeOutCode(x1, y1);
}
}
}
if (accept) {
// printf( "x0 = %d, y0 = %d, x1 = %d, y1 =%d \n",x0, y0, x1, y1);
dDrawLineSegment(x0, y0, x1, y1, iColor, A);
}
}
// -----------------------------------------------------------------------
int DrawClippedLineSegment(complex double z0, complex double z1, unsigned char iColor, unsigned char A[]){
double x0 = creal(z0);
double y0 = cimag(z0);
double x1 = creal(z1);
double y1 = cimag(z1);
dDrawLine( x0, y0, x1, y1, iColor, A);
return 0;
}
//------------------mapping between world and integer coordinate -----------------------------------------------------
// 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;
}
// ****************** DYNAMICS = trap tests ( target sets) ****************************
// 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
//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;
//}
// compute alfa fixed point
// https://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings#Period-1_points_(fixed_points)
complex double GiveAlfa(complex double c)
{
// d=1-4c
// alfa = (1-sqrt(d))/2
return (1.0-csqrt(1.0 - 4.0*c))/2.0 ;
}
// *************************************************************************************************
// ********************************* critical orbit ************************************************
// *********************************************************************************************
// fill array
// uses global var : ...
int DrawImage_CriticalOrbit (unsigned char A[])
{
fprintf(stderr, "draw critical orbit \n");
int i = 0; // iteration = number of the point
int iMax = child_period*iterMax_cr_orbit;
complex double z = 0.0;
for (i = 0; i < iMax; ++i){
fprintf(stderr, "\t%d / %d\r", i, iMax); // info
ColorPixel_d(z, 255, A ); // draw point and check if point is outside image
z = z*z+c; // forward iteration : complex quadratic polynomial
}
return 0;
}
// fill array
// uses global var : ...
int DrawImage_CriticalOrbit_Approx (unsigned char A[])
{
fprintf(stderr, "draw critical orbit \n");
int i = 0; // iteration = number of the point
int iMax = child_period*iterMax_cr_orbit;
complex double z = 0.0;
for (i = 0; i < iMax; ++i){
fprintf(stderr, "\t%d / %d\r", i, iMax); // info
ColorPixel_d(z, 255, A ); // draw point and check if point is outside image
z = z*z+c; // forward iteration : complex quadratic polynomial
}
// aproximation and info
double t = ((double)e_angle_numerator0)/e_angle_denominator; // first external angle in turns
for (i = 0; i < child_period; ++i){
printf (" i = %d \t t = %.16f \tz = %.16f %+.16f*I \t distance to fixed point = %.16f = %.16f*PixelWidth = %.16f*ImageWidth\n ", i, t, creal(z), cimag(z), cabs(a - z), cabs(a-z)/PixelWidth, cabs(a-z)/ImageWidth );
// approximate when landing point is known
DrawClippedLineSegment(z, a, 155, A) ;
// next angle
t *= 2.0; // t = 2*t angle doubling map
if (t > 1.0) {t--;} // t = t modulo 1
// next z
z = z*z+c; // forward iteration : complex quadratic polynomial
}
return 0;
}
// *************************************************************************************************
// ********************************* external ray ************************************************
// *********************************************************************************************
// check if period p of n/d is proper under period doubling map
// p times (2*n) mod d should give n
int TestPeriod(int n, int d, int period){
int p;
int pMax = period;
int nn = n; // termporary value
if (n<=0 || d<=0 || period <= 0) {
fprintf(stderr, " input should be positive \n");
return 1;
}
for(p=0; p<pMax; p++)
nn *=2;
nn = nn % d; // remainder
if (nn-n ==0)
return 1; // true
return 0; // false
}
// dot product
double dot(complex double a, complex double b){
// ax*bx + ay*by
return creal(a)*creal(b) + cimag(a)*cimag(b);
}
int InverseIterationAndDrawSegment (int per, complex double zz[per], unsigned char A[]){
int p;
int pMax = per;
complex double z;
complex double zn; // next vsalue : zn = f^(-1)(z)
// compute preimages
for(p=0; p<pMax; ++p){
z = zz[(p+1) % per]; //read point from the ray p+1
zn = csqrt(z-c); // inverse iteration = preimage; fc(z) maps z_{l,j} to z_{l-1,j+1}} so f^{-1} map to z_{l+1, j-1}
// choose correct preimage using dot product
if (dot(zz[p],zn) < 0 ) zn = - zn;
// draw segment of ray p
DrawClippedLineSegment(zz[p],zn, 255, A);
zz[p] = zn; //save
}
return 0;
}
// uses global var : ...
/*
"In the dynamic plane, external rays can be drawn by backwards iteration. It is most effective for a periodic or preperiodic angle.
You must keep track of points on the finite collection of rays with angles phi ,2phi ,4phi ...
Say z_{l,j} corresponds to a radius R^{1/(2l)} and the angle 2j*phi .
Then fc(z) maps z_{l,j} to z_{l-1,j+1}}
This point, which was constructed before, has two preimages under {\displaystyle f_{c}(z)} {\displaystyle f_{c}(z)} .
The one that is closer to {\displaystyle z_{l-1,j}} {\displaystyle z_{l-1,j}} is the correct one.
This criterion was proved by Thierry Bousch. The ray will look better when you introduce intermediate points." Wolf Jung
this procedure works only for periodic angles
*/
int DrawExternalDynamicRaysBI (int n, int m, int period, int iMax, unsigned char A[])
{ // In the dynamic plane, external rays can be drawn by backwards iteration. This procedure is only for the periodic angle
// https://commons.wikimedia.org/wiki/File:Backward_Iteration.svg
fprintf(stderr, "draw external ray \n");
int i = 0; // iteration = number of points
//int iMax = 10000000;
double r = 10000.0; // very big radius = near infinity where z=w so one can swith to dynamical plane ( Boettcher conjugation )
complex double zz[period]; // zz is an array of z
double t;
int p;
int pMax = period; // number of rays to draw
int result ;
// check if period is proper
// p times (2*n) mod m should give n
if (TestPeriod(n, m, period)!=1 ) {
fprintf(stderr, "bad input to the TestPeriod from external ray : error \n");
return 1;
}
fprintf(stderr, "good input to the TestPeriod from external ray \n");
// to do
// find landing point : periodic point
// compute how many times ray rotates around it's landing point
// compute distance to tha landing point of last point of the ray ( along the ray ? = infinity so to the circle with fixed radius around landing point)
// initial points on rays
t = (double)n/m; // first external angle in turns
for(p=0; p<pMax; ++p){
// initial point
zz[p] = r*cexp(2.0*I * M_PI * t ); // Euler's formula
// next angle
t *= 2.0; // t = 2*t angle doubling map
if (t > 1.0) t--; // t = t modulo 1
}
for (i = 0; i < iMax; ++i){
// inverse iteration of complex quadratic polynomial: z = csqrt(z-c) with proper choose of preimage for one point on every ray
fprintf(stderr, "\t%d / %d\r", i, iMax); // info
result = InverseIterationAndDrawSegment(period, zz, A);
if (result != 0) {fprintf (stderr, " lost precision\n"); return 1;}
}
fprintf(stderr, "end of drawing external ray \n");
return 0;
}
// uses global var : ...
/*
"In the dynamic plane, external rays can be drawn by backwards iteration. It is most effective for a periodic or preperiodic angle.
You must keep track of points on the finite collection of rays with angles phi ,2phi ,4phi ...
Say z_{l,j} corresponds to a radius R^{1/(2l)} and the angle 2j*phi .
Then fc(z) maps z_{l,j} to z_{l-1,j+1}}
This point, which was constructed before, has two preimages under {\displaystyle f_{c}(z)} {\displaystyle f_{c}(z)} .
The one that is closer to {\displaystyle z_{l-1,j}} {\displaystyle z_{l-1,j}} is the correct one.
This criterion was proved by Thierry Bousch. The ray will look better when you introduce intermediate points." Wolf Jung
this procedure works only for periodic angles
*/
int DrawExternalDynamicRaysBI_approx (int n, int m, int period, int iMax, unsigned char A[])
{ // In the dynamic plane, external rays can be drawn by backwards iteration. This procedure is only for the periodic angle
// https://commons.wikimedia.org/wiki/File:Backward_Iteration.svg
fprintf(stderr, "draw external ray \n");
int i = 0; // iteration = number of points
//int iMax = 10000000;
double r = 10000.0; // very big radius = near infinity where z=w so one can swith to dynamical plane ( Boettcher conjugation )
complex double zz[period]; // zz is an array of z
double t;
int p;
int pMax = period; // number of rays to draw
int result ;
// check if period is proper
// p times (2*n) mod m should give n
if (TestPeriod(n, m, period)!=1 ) {
fprintf(stderr, "bad input to the TestPeriod from external ray : error \n");
return 1;
}
fprintf(stderr, "good input to the TestPeriod from external ray \n");
// to do
// find landing point : periodic point
// compute how many times ray rotates around it's landing point
// compute distance to tha landing point of last point of the ray ( along the ray ? = infinity so to the circle with fixed radius around landing point)
// initial points on rays
t = (double)n/m; // first external angle in turns
for(p=0; p<pMax; ++p){
// initial point
zz[p] = r*cexp(2.0*I * M_PI * t ); // Euler's formula
// next angle
t *= 2.0; // t = 2*t angle doubling map
if (t > 1.0) t--; // t = t modulo 1
}
for (i = 0; i < iMax; ++i){
// inverse iteration of complex quadratic polynomial: z = csqrt(z-c) with proper choose of preimage for one point on every ray
fprintf(stderr, "\t%d / %d\r", i, iMax); // info
result = InverseIterationAndDrawSegment(period, zz, A);
if (result != 0) {fprintf (stderr, " lost precision\n"); return 1;}
}
// print and approximate if one knows landing point !
printf("last point of the ray\n");
t = (double)n/m; // first external angle in turns
for(p=0; p<pMax; ++p){
printf (" p = %d \t t = %.16f \tz = %.16f %+.16f*I \t distance to landing point = %.16f = %.16f* PixelWidth = %.16f*ImageWidth\n ", p, t, creal(zz[p]), cimag(zz[p]), cabs(a - zz[p]), cabs(a- zz[p])/PixelWidth, cabs(a- zz[p])/ImageWidth );
// next angle
t *= 2.0; // t = 2*t angle doubling map
if (t > 1.0) {t--;} // t = t modulo 1
// approximation when landing point is known
DrawClippedLineSegment(zz[p], a, 155, A) ;
}
fprintf(stderr, "end of drawing external ray \n");
return 0;
}
// *******************************************************************************************
// ********************************** save A array to pgm file ****************************
// *********************************************************************************************
int SaveArray2PGMFile( unsigned char A[], double k, 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 [200]; /* name of file */
snprintf(name, sizeof name, "z%.0f", k); /* */
char *filename =strcat(name,".pgm");
// 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(A,iSize,1,fp); // write array with image data bytes to the file in one step
fclose(fp);
// info
fprintf(stderr, "File %s saved ", filename);
if (comment == NULL || strlen(comment) ==0)
fprintf(stderr,"\n");
else fprintf (stderr, ". Comment = %s \n", comment);
return 0;
}
int PrintInfoAboutProgam()
{
// display info messages
printf ("\n\nNumerical 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 ) \n", creal(c), cimag(c));
printf ("fixed point alfa z = a = ( %.16f ; %.16f ) \n", creal(a), cimag(a));
printf ("iSize = %d\n", iSize);
printf ("Image Width = %.16f in world coordinate\n", ZxMax - ZxMin);
// https://en.wikipedia.org/wiki/Aspect_ratio_(image)
printf ("Image aspect ratio = %.16f\n", ratio);
printf ("PixelWidth = %.16f \n", PixelWidth);
printf ("\nexternal ray \n");
printf ("Maximal number of iterations = iterMax_ray = %ld \n", iterMax_ray);
printf ("Maximal number of iterations = iterMax_cr_orbit = %ld \n", iterMax_cr_orbit);
printf ("MinDistanceToLandingPoint = %.16f \n", MinDistanceToLandingPoint);
printf ("\n plane : \n");
printf (" radius = %.16f \n", plane_radius);
printf (" center = z = %.16f %+.16f*I )\n ", creal(plane_center), cimag(plane_center) );
// image corners in world coordinate
// center and radius
// center and zoom
// GradientRepetition
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 diplayed in the console
return 0;
}
// *****************************************************************************
//;;;;;;;;;;;;;;;;;;;;;; setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
// **************************************************************************************
int setup ()
{
fprintf (stderr, "setup start ");
c = -0.690059870015044 +0.276026482784614*I; // t = 5/11
a = GiveAlfa(c);
t_denominator = child_period;
plane_center = a;
plane_radius = 0.9;
/* 2D array ranges */
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].
//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].
ratio = (double) iWidth / iHeight;
SetPlane(plane_center , plane_radius, ratio); //
/* Pixel sizes */
PixelWidth = (ZxMax - ZxMin) / ixMax; // ixMax = (iWidth-1) step between pixels in world coordinate
PixelHeight = (ZyMax - ZyMin) / iyMax;
//((ZxMax - ZxMin) / (ZyMax - ZyMin)) / ((double) iWidth / (double) iHeight); // it should be 1.000 ...
//ER2 = ER * ER; // for numerical optimisation in iteration
//lnER = log(EscapeRadius); // ln(ER)
/* 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 = 6.0*iWidth/2000.0 ; // measured in pixels ( when iWidth = 2000) ; such function is stable when iWidth is changing
//distanceMax = BoundaryWidth*PixelWidth; // distance to the boundary from exterior
//InternalSiegelDiscRadius = GiveInternalSiegelDiscRadius(c,a);
//ExternalSiegelDiscRadius = GiveExternalSiegelDiscRadius(c,a);
fprintf (stderr, "and end\n");
return 0;
} // ;;;;;;;;;;;;;;;;;;;;;;;;; end of the setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
int end(){
// printf (" 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 ();
DrawImage_CriticalOrbit_Approx(data);
SaveArray2PGMFile (data, child_period*iterMax_cr_orbit, "critical orbit");
DrawExternalDynamicRaysBI_approx(e_angle_numerator0,e_angle_denominator,child_period, child_period*iterMax_ray, data);
SaveArray2PGMFile (data, child_period*iterMax_ray, "ray");
end();
return 0;
}
bash source code
#!/bin/bash
# script file for BASH
# which bash
# save this file as j.sh
# chmod +x j.sh
# ./j.sh
# checked in https://www.shellcheck.net/
printf "make pgm files \n"
gcc e.c -lm -Wall -march=native
if [ $? -ne 0 ]
then
echo ERROR: compilation failed !!!!!!
exit 1
fi
printf "run the compiled program\n"
time ./a.out > e.txt
export OMP_DISPLAY_ENV="FALSE"
printf "change Image Magic settings\n"
export MAGICK_WIDTH_LIMIT=100MP
export MAGICK_HEIGHT_LIMIT=100MP
printf "convert all pgm files to png using Image Magic v 6 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
convert "${b}".pgm "${b}".png
echo "$file"
done
printf "delete all pgm files \n"
rm ./*.pgm
echo OK
printf "info about software \n"
bash --version
make -v
gcc --version
convert -version
convert -list resource
# end
make
all:
chmod +x d.sh
./d.sh
Tu run the program simply
make
text output
Critical orbit: last point i = 0 t = 0.1665852467024914 z = -0.5304160097901404 +0.0994122341565338*I distance to fixed point = 0.0654663145836885 = 72.7039793626629347*PixelWidth = 0.0363701747687158*ImageWidth i = 1 t = 0.3331704934049829 z = -0.4186015188733432 +0.1705668016533505*I distance to fixed point = 0.0679766738222065 = 75.4918727614393248*PixelWidth = 0.0377648187901147*ImageWidth i = 2 t = 0.6663409868099658 z = -0.5439256722382275 +0.1332274383016925*I distance to fixed point = 0.0646321879634904 = 71.7776354105652103*PixelWidth = 0.0359067710908280*ImageWidth i = 3 t = 0.3326819736199316 z = -0.4119542834116676 +0.1310948349069639*I distance to fixed point = 0.0684928021734420 = 76.0650619692835903*PixelWidth = 0.0380515567630233*ImageWidth i = 4 t = 0.6653639472398631 z = -0.5375393941331076 +0.1680163252384757*I distance to fixed point = 0.0638525268049897 = 70.9117783795413743*PixelWidth = 0.0354736260027721*ImageWidth i = 5 t = 0.3307278944797263 z = -0.4293407553166969 +0.0953956954382913*I distance to fixed point = 0.0678845467235295 = 75.3895605001863629*PixelWidth = 0.0377136370686275*ImageWidth i = 6 t = 0.6614557889594526 z = -0.5148267245472875 +0.1941119629177389*I distance to fixed point = 0.0637630457040305 = 70.8124046457539009*PixelWidth = 0.0354239142800170*ImageWidth i = 7 t = 0.3229115779189051 z = -0.4626927678547330 +0.0761584306558459*I distance to fixed point = 0.0669173736220485 = 74.3154610391527228*PixelWidth = 0.0371763186789158*ImageWidth i = 8 t = 0.6458231558378102 z = -0.4817753791499315 +0.2055505726333618*I distance to fixed point = 0.0647161055837454 = 71.8708305899483975*PixelWidth = 0.0359533919909697*ImageWidth i = 9 t = 0.2916463116756205 z = -0.5002033919698867 +0.0779680726547672*I distance to fixed point = 0.0661412825503758 = 73.4535687878895942*PixelWidth = 0.0367451569724310*ImageWidth i = 10 t = 0.5832926233512410 z = -0.4459354570303629 +0.1980266939700758*I distance to fixed point = 0.0664115866423103 = 73.7537564988768537*PixelWidth = 0.0368953259123946*ImageWidth last point of the external ray p = 0 t = 0.1665852467024914 z = -0.4170188752901364 +0.1506018053351233*I distance to landing point = 0.0634786084515321 = 70.4965212747848113* PixelWidth = 0.0352658935841845*ImageWidth p = 1 t = 0.3331704934049829 z = -0.5388360314336247 +0.1504188918297714*I distance to landing point = 0.0598567181450451 = 66.4742108733029085* PixelWidth = 0.0332537323028029*ImageWidth p = 2 t = 0.6663409868099658 z = -0.4223414442666910 +0.1139242453319452*I distance to landing point = 0.0634130274946685 = 70.4236899788012778* PixelWidth = 0.0352294597192603*ImageWidth p = 3 t = 0.3326819736199316 z = -0.5246663081382248 +0.1797966221587810*I distance to landing point = 0.0594420895659870 = 66.0137428013377558* PixelWidth = 0.0330233830922150*ImageWidth p = 4 t = 0.6653639472398631 z = -0.4471119604606156 +0.0873600228620700*I distance to landing point = 0.0626732135268017 = 69.6020854667092408* PixelWidth = 0.0348184519593343*ImageWidth p = 5 t = 0.3307278944797263 z = -0.4977825384230178 +0.1979070606045213*I distance to landing point = 0.0598243260724863 = 66.4382376771666827* PixelWidth = 0.0332357367069368*ImageWidth p = 6 t = 0.6614557889594526 z = -0.4814396190918361 +0.0789971247905459*I distance to landing point = 0.0618923167109339 = 68.7348561695315681* PixelWidth = 0.0343846203949633*ImageWidth p = 7 t = 0.3229115779189051 z = -0.4645163089093985 +0.1999617914428644*I distance to landing point = 0.0610265350331413 = 67.7733575173608500* PixelWidth = 0.0339036305739674*ImageWidth p = 8 t = 0.6458231558378102 z = -0.5142691868062246 +0.0902554562208356*I distance to landing point = 0.0612639546468155 = 68.0370251883245771* PixelWidth = 0.0340355303593420*ImageWidth p = 9 t = 0.2916463116756205 z = -0.4337331208976685 +0.1831952826301759*I distance to landing point = 0.0625217917185160 = 69.4339231362852871* PixelWidth = 0.0347343287325089*ImageWidth p = 10 t = 0.5832926233512410 z = -0.5354959614261718 +0.1171107594494743*I distance to landing point = 0.0605997409390414 = 67.2993789650798533* PixelWidth = 0.0336665227439119*ImageWidth Numerical approximation of Julia set for fc(z)= z^2 + c parameter c = ( -0.6900598700150440 ; 0.2760264827846140 ) fixed point alfa z = a = ( -0.4797464868072486 ; 0.1408662784207147 ) iSize = 4000000 Image Width = 1.8000000000000000 in world coordinate Image aspect ratio = 1.0000000000000000 PixelWidth = 0.0009004502251126 external ray Maximal number of iterations = iterMax_ray = 2000000000000 Maximal number of iterations = iterMax_cr_orbit = 1000000000000 MinDistanceToLandingPoint = 1.0000000000000000 plane : radius = 0.9000000000000000 center = z = -0.4797464868072486 +0.1408662784207147*I ) ratio of image = 1.000000 ; it should be 1.000 ... gcc version: 11.3.0 info about software GNU bash, version 5.1.16(1)-release (x86_64-pc-linux-gnu) Copyright (C) 2020 Free Software Foundation, Inc. License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html> This is free software; you are free to change and redistribute it. There is NO WARRANTY, to the extent permitted by law. GNU Make 4.3 Built for x86_64-pc-linux-gnu Copyright (C) 1988-2020 Free Software Foundation, Inc. License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html> This is free software: you are free to change and redistribute it. There is NO WARRANTY, to the extent permitted by law. gcc (Ubuntu 11.3.0-1ubuntu1~22.04) 11.3.0 Copyright (C) 2021 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. Version: ImageMagick 6.9.11-60 Q16 x86_64 2021-01-25 https://imagemagick.org Copyright: (C) 1999-2021 ImageMagick Studio LLC License: https://imagemagick.org/script/license.php Features: Cipher DPC Modules OpenMP(4.5) Delegates (built-in): bzlib djvu fftw fontconfig freetype heic jbig jng jp2 jpeg lcms lqr ltdl lzma openexr pangocairo png tiff webp wmf x xml zlib Resource limits: Width: 1MP Height: 1MP List length: unlimited Area: 128MP Memory: 256MiB Map: 512MiB Disk: 10GiB File: 768 Thread: 8 Throttle: 0 Time: unlimited
references
Items portrayed in this file
depicts
some value
15 March 2023
File history
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 13:27, 15 March 2023 | 2,000 × 2,000 (27 KB) | Soul windsurfer | Uploaded own work with UploadWizard |
File usage
No pages on the English Wikipedia use this file (pages on other projects are not listed).
Global file usage
The following other wikis use this file:
Metadata
This file contains additional information, probably added from the digital camera or scanner used to create or digitize it.
If the file has been modified from its original state, some details may not fully reflect the modified file.
PNG file comment |
|
---|---|
File change date and time | 12:23, 15 March 2023 |