קובץ:Parabolic Julia set for internal angle 1 over 7 with flower trap and critical orbit.png
תוכן הדף אינו נתמך בשפות אחרות.
מראה
מתוך ויקיפדיה, האנציקלופדיה החופשית
גודל התצוגה המקדימה הזאת: 600 × 600 פיקסלים. רזולוציות אחרות: 240 × 240 פיקסלים | 480 × 480 פיקסלים | 1,000 × 1,000 פיקסלים.
לקובץ המקורי (1,000 × 1,000 פיקסלים, גודל הקובץ: 100 ק"ב, סוג MIME: image/png)
זהו קובץ שמקורו במיזם ויקישיתוף. תיאורו בדף תיאור הקובץ המקורי (בעברית) מוצג למטה. |
תוכן עניינים
תקציר
תיאורParabolic Julia set for internal angle 1 over 7 with flower trap and critical orbit.png |
English: Numerical approximation of parabolic Julia set for :
|
תאריך יצירה | |
מקור | Own program which uses the algorithm and code by Wolf Jung |
יוצר | Adam majewski |
רישיון
אני, בעל זכויות היוצרים על עבודה זו, מפרסם בזאת את העבודה תחת הרישיון הבא:
הקובץ הזה מתפרסם לפי תנאי רישיון קריאייטיב קומונז ייחוס-שיתוף זהה 3.0 לא מותאם.
- הנכם רשאים:
- לשתף – להעתיק, להפיץ ולהעביר את העבודה
- לערבב בין עבודות – להתאים את העבודה
- תחת התנאים הבאים:
- ייחוס – יש לתת ייחוס הולם, לתת קישור לרישיון, ולציין אם נעשו שינויים. אפשר לעשות את זה בכל צורה סבירה, אבל לא בשום צורה שמשתמע ממנה שמעניק הרישיון תומך בך או בשימוש שלך.
- שיתוף זהה – אם תיצרו רמיקס, תשנו, או תבנו על החומר, חובה עליכם להפיץ את התרומות שלך לפי תנאי רישיון זהה או תואם למקור.
C src code
Src code was formatted with Emacs
/*
Adam Majewski
fraktal.republika.pl
c console progam using
* symmetry
* openMP
draw parabolic Julia set
and saves it to pgm file
gcc t.c -lm -Wall -fopenmp -march=native
time ./a.out
*/
#include <stdio.h>
#include <stdlib.h> // malloc
#include <string.h> // strcat
#include <math.h> // M_PI; needs -lm also
#include <complex.h>
#include <omp.h> // OpenMP; needs also -fopenmp
/* --------------------------------- global variables and constans ------------------------------------------------------------ */
//unsigned int period = 7; // period of secondary component joined by root point
#define period 7
// unsigned int denominator; denominator = period;
double InternalAngle;
unsigned char Colors[period]; //={255,230,180, 160,140,120,100}; // NumberOfPetal of colors >= period
unsigned char iExterior = 245;
unsigned char iPetal = 255;
// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1
unsigned int ix, iy; // var
unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
unsigned int ixMax ; //
unsigned int iWidth ; // horizontal dimension of array
unsigned int ixAxisOfSymmetry ; //
unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
unsigned int iyMax ; //
unsigned int iyAxisOfSymmetry ; //
unsigned int iyAbove ; // var, measured from 1 to (iyAboveAxisLength -1)
unsigned int iyAboveMin = 1 ; //
unsigned int iyAboveMax ; //
unsigned int iyAboveAxisLength ; //
unsigned int iyBelowAxisLength ; //
unsigned int iHeight = 2000; // odd NumberOfPetal !!!!!! = (iyMax -iyMin + 1) = iyAboveAxisLength + iyBelowAxisLength +1
// The size of array has to be a positive constant integer
unsigned int iSize ; // = iWidth*iHeight;
// memmory 1D arrays
unsigned char *data;
unsigned char *edge;
// unsigned int i; // var = index of 1D array
unsigned int iMin = 0; // Indexes of array starts from 0 not 1
unsigned int iMax ; // = i2Dsize-1 =
// The size of array has to be a positive constant integer
// unsigned int i1Dsize ; // = i2Dsize = (iMax -iMin + 1) = ; 1D array with the same size as 2D array
/* world ( double) coordinate = dynamic plane */
const double ZxMin=-1.5;
const double ZxMax=1.5;
const double ZyMin=-1.5;
const double ZyMax=1.5;
double PixelWidth; // =(ZxMax-ZxMin)/iXmax;
double PixelHeight; // =(ZyMax-ZyMin)/iYmax;
double ratio ;
// complex numbers of parametr plane
double Cx; // c =Cx +Cy * i
double Cy;
double complex c; //
double complex alfa; // alfa fixed point alfa=f(alfa)
double alfax,alfay;
unsigned long int iterMax = 10000; //iHeight*100;
// target set for escaping points is a exterior of circle with center in origin
double ER = 2.0; // Escape Radius for bailout test
double ER2;
/* trap triangle parameters */
double d1=-0.4; /* length */
double d2=0.6; /* width */
/* ------------------------------------------ functions -------------------------------------------------------------*/
/* find c in component of Mandelbrot set
uses code by Wolf Jung from program Mandel
see function mndlbrot::bifurcate from mandelbrot.cpp
http://www.mndynamics.com/indexp.html
*/
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int iPeriod)
{
//0 <= InternalRay<= 1
//0 <= InternalAngleInTurns <=1
double t = InternalAngleInTurns *2*M_PI; // from turns to radians
double R2 = InternalRadius * InternalRadius;
//double Cx, Cy; /* C = Cx+Cy*i */
switch ( iPeriod ) // of component
{
case 1: // main cardioid
Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4;
Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4;
break;
case 2: // only one component
Cx = InternalRadius * 0.25*cos(t) - 1.0;
Cy = InternalRadius * 0.25*sin(t);
break;
// for each period there are 2^(period-1) roots.
default: // higher periods : to do
Cx = 0.0;
Cy = 0.0;
break; }
return Cx + Cy*I;
}
/*
http://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings
z^2 + c = z
z^2 - z + c = 0
ax^2 +bx + c =0 // general form of quadratic equation
so :
a=1
b =-1
c = c
so :
The discriminant is the d=b^2- 4ac
d=1-4c = dx+dy*i
r(d)=sqrt(dx^2 + dy^2)
sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx +- sy*i
x1=(1+sqrt(d))/2 = beta = (1+sx+sy*i)/2
x2=(1-sqrt(d))/2 = alfa = (1-sx -sy*i)/2
alfa : attracting when c is in main cardioid of Mandelbrot set, then it is in interior of Filled-in Julia set,
it means belongs to Fatou set ( strictly to basin of attraction of finite fixed point )
*/
// uses global variables :
// ax, ay (output = alfa(c))
double complex GiveAlfaFixedPoint(double complex c)
{
double dx, dy; //The discriminant is the d=b^2- 4ac = dx+dy*i
double r; // r(d)=sqrt(dx^2 + dy^2)
double sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i
double ax, ay;
// d=1-4c = dx+dy*i
dx = 1 - 4*creal(c);
dy = -4 * cimag(c);
// r(d)=sqrt(dx^2 + dy^2)
r = sqrt(dx*dx + dy*dy);
//sqrt(d) = s =sx +sy*i
sx = sqrt((r+dx)/2);
sy = sqrt((r-dx)/2);
// alfa = ax +ay*i = (1-sqrt(d))/2 = (1-sx + sy*i)/2
ax = 0.5 - sx/2.0;
ay = sy/2.0;
return ax+ay*I;
}
// colors of components interior = shades of gray
int InitColors()
{
int i;
int iMax = period; // uses global var period and Colors
unsigned int iStep;
iStep=150/period;
for (i = 1; i <= iMax; ++i)
{Colors[i-1] = iExterior -i*iStep;
printf("i= %d color = %i \n",i-1, Colors[i-1]);}
return 0;
}
/* -------------------------------------------------- SETUP --------------------------------------- */
int setup()
{
/* 2D array ranges */
if (!(iHeight % 2)) iHeight+=1; // it sholud be even NumberOfPetal (variable % 2) or (variable & 1)
iWidth = iHeight;
iSize = iWidth*iHeight; // size = NumberOfPetal of points in array
// iy
iyMax = iHeight - 1 ; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].
iyAboveAxisLength = (iHeight -1)/2;
iyAboveMax = iyAboveAxisLength ;
iyBelowAxisLength = iyAboveAxisLength; // the same
iyAxisOfSymmetry = iyMin + iyBelowAxisLength ;
// ix
ixMax = iWidth - 1;
/* 1D array ranges */
// i1Dsize = i2Dsize; // 1D array with the same size as 2D array
iMax = iSize-1; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].
InternalAngle = 1.0/((double) period); // angle for 3/7 is different then angle for 2/7
c = GiveC(InternalAngle, 1.0, 1) ;
alfa = GiveAlfaFixedPoint(c);
Cx=creal(c);
Cy=cimag(c);
alfax=creal(alfa);
alfay=cimag(alfa);
/* Pixel sizes */
PixelWidth = (ZxMax-ZxMin)/ixMax; // ixMax = (iWidth-1) step between pixels in world coordinate
PixelHeight = (ZyMax-ZyMin)/iyMax;
ratio = ((ZxMax-ZxMin)/(ZyMax-ZyMin))/((float)iWidth/(float)iHeight); // it should be 1.000 ...
// for numerical optimisation
ER2 = ER * ER;
/* create dynamic 1D arrays for colors ( shades of gray ) */
data = malloc( iSize * sizeof(unsigned char) );
edge = malloc( iSize * sizeof(unsigned char) );
if (edge == NULL || edge == NULL )
{
fprintf(stderr," Could not allocate memory\n");
return 1;
}
else fprintf(stderr," memory is OK \n");
InitColors();
return 0;
}
// from screen to world coordinate ; linear mapping
// uses global cons
double GiveZx(unsigned int ix)
{ return (ZxMin + ix*PixelWidth );}
// uses global cons
double GiveZy(unsigned int iy)
{ return (ZyMax - iy*PixelHeight);} // reverse y axis
/*
Checks to which petal ( = part of target set) point of interior falls
bail-out test = check if points escapes ) return 0
if not escapes find what iteration it need to fall into trap ( target set inside one petal)
*/
/*
uses code by Wolf Jung from program Mandel
see function mndlbrot::pixcolor from ....cpp
http://www.mndynamics.com/indexp.html
uint mndlbrot::pixcolor(double x, double y)
{ uint j; double a = A, b = b, x1, y1, u, v;
for (j = 1; j <= 1000*maxiter; j++) //larger number !!!!!!!!!!
{ u = x*x; v = y*y;
if (u + v <= bailout) { y = 2*x*y + b; x = u - v + a; }
else break;
x1 = x - 0.311744900929367; y1 = y -0.390915741234015;
if (x1 + y1 > -0.5 && y1 > x1 && y1 < 0.4*x1)
return 9 + j % 7;
}
return 0;
}
*/
unsigned char GiveColor(unsigned int ix, unsigned int iy)
{
double Zx, Zy; // Z= Zx+ZY*i;
double Zx2, Zy2;
double xt, yt;
int j;
// from screen to world coordinate
Zx = GiveZx(ix);
Zy = GiveZy(iy);
for (j = 1; j <= iterMax; j++) //larg number of iteration s
{ Zx2 = Zx*Zx;
Zy2 = Zy*Zy;
// bailout test
if (Zx2 + Zy2 > ER2) return iExterior; // if escaping stop iteration
// if not escaping iterate
// Z(n+1) = Zn * Zn + C
Zy = 2*Zx*Zy + Cy;
Zx = Zx2 - Zy2 + Cx;
// translation near fixed point alfa
xt = Zx - alfax;
yt = Zy -alfay;
// check if point z is inside triangle around arm of critical orbit
if (xt + yt > d1 && yt > xt && yt < d2*xt) return Colors[ j % period]; //
}
return iPetal; // not escaping and not in attracting target set
}
/* ----------- array functions -------------- */
/* gives position of 2D point (iX,iY) in 1D array ; uses also global variable iWidth */
unsigned int Give_i(unsigned int ix, unsigned int iy)
{ return ix + iy*iWidth; }
// ix = i % iWidth;
// iy = (i- ix) / iWidth;
// i = Give_i(ix, iy);
// plots raster point (ix,iy)
int PlotPoint(unsigned int ix, unsigned int iy, unsigned char iColor)
{
unsigned i; /* index of 1D array */
i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
data[i] = iColor;
return 0;
}
// fill array
// uses global var : ...
// scanning complex plane
int FillArray(unsigned char data[] )
{
unsigned int ix, iy; // pixel coordinate
// for all pixels of image
for(iy = iyMin; iy<=iyMax; ++iy)
{ printf(" %d z %d\n", iy, iyMax); //info
for(ix= ixMin; ix<=ixMax; ++ix) PlotPoint(ix, iy, GiveColor(ix, iy) ); //
}
return 0;
}
// fill array using symmetry of image
// uses global var : ...
int FillArraySymmetric(unsigned char data[] )
{
unsigned char Color; // gray from 0 to 255
printf("axis of symmetry \n");
iy = iyAxisOfSymmetry;
#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)
for(ix=ixMin;ix<=ixMax;++ix) {//printf(" %d from %d\n", ix, ixMax); //info
PlotPoint(ix, iy, GiveColor(ix, iy));
}
/*
The use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.
The use of ‘private(variable, variable2)’ specifies that these variables should have a seperate instance in each thread.
*/
#pragma omp parallel for schedule(dynamic) private(iyAbove,ix,iy,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)
// above and below axis
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove)
{printf(" %d from %d\r", iyAbove, iyAboveMax); //info
for(ix=ixMin; ix<=ixMax; ++ix)
{ // above axis compute color and save it to the array
iy = iyAxisOfSymmetry + iyAbove;
Color = GiveColor(ix, iy);
PlotPoint(ix, iy, Color );
// below the axis only copy Color the same as above without computing it
PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color );
}
}
return 0;
}
int AddBoundaries(unsigned char data[])
{
unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
unsigned int i; /* index of 1D array */
/* sobel filter */
unsigned char G, Gh, Gv;
printf(" find boundaries in data array using Sobel filter\n");
#pragma omp parallel for schedule(dynamic) private(i,iY,iX,Gv,Gh,G) shared(iyMax,ixMax, ER2)
for(iY=1;iY<iyMax-1;++iY){
for(iX=1;iX<ixMax-1;++iX){
Gv= data[Give_i(iX-1,iY+1)] + 2*data[Give_i(iX,iY+1)] + data[Give_i(iX-1,iY+1)] - data[Give_i(iX-1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[Give_i(iX+1,iY-1)];
Gh= data[Give_i(iX+1,iY+1)] + 2*data[Give_i(iX+1,iY)] + data[Give_i(iX-1,iY-1)] - data[Give_i(iX+1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[Give_i(iX-1,iY-1)];
G = sqrt(Gh*Gh + Gv*Gv);
i= Give_i(iX,iY); /* compute index of 1D array from indices of 2D array */
if (G==0) {edge[i]=255;} /* background */
else {edge[i]=0;} /* boundary */
}
}
// copy boundaries from edge array to data array
//for(iY=1;iY<iyMax-1;++iY){
// for(iX=1;iX<ixMax-1;++iX){i= Give_i(iX,iY); if (edge[i]==0) data[i]=0;}}
return 0;
}
// Check Orientation of image : first quadrant in upper right position
// uses global var : ...
int CheckOrientation(unsigned char data[] )
{
unsigned int ix, iy; // pixel coordinate
double Zx, Zy; // Z= Zx+ZY*i;
unsigned i; /* index of 1D array */
for(iy=iyMin;iy<=iyMax;++iy)
{
Zy = GiveZy(iy);
for(ix=ixMin;ix<=ixMax;++ix)
{
// from screen to world coordinate
Zx = GiveZx(ix);
i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
if (Zx>0 && Zy>0) data[i]=255-data[i]; // check the orientation of Z-plane by marking first quadrant */
}
}
return 0;
}
int CopyBoundaries()
{
unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
unsigned int i; /* index of 1D array */
printf("copy boundaries from edge array to data array \n");
for(iY=1;iY<iyMax-1;++iY)
for(iX=1;iX<ixMax-1;++iX)
{i= Give_i(iX,iY); if (edge[i]==0) data[i]=0;}
return 0;
}
int MarkTrap()
{
unsigned int ix, iy; // pixel coordinate
double Zx, Zy; // Z= Zx+ZY*i;
double xt, yt;
unsigned i; /* index of 1D array */
for(iy=iyMin;iy<=iyMax;++iy)
{
Zy = GiveZy(iy);
for(ix=ixMin;ix<=ixMax;++ix)
{
// from screen to world coordinate
Zx = GiveZx(ix);
i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
// translation near fixed point alfa
xt = Zx - alfax;
yt = Zy -alfay;
// check if point z is inside triangle around arm of critical orbit
if (xt + yt > d1 && yt > xt && yt < d2*xt) data[i]=255-data[i]; // mark the trap
}
}
return 0;
}
int DrawCriticalOrbit(unsigned int IterMax)
{
unsigned int ix, iy; // pixel coordinate
double Zx=0.0;
double Zy=0.0; // Z= Zx+ZY*i;
double Zx2=0.0;
double Zy2=0.0;
unsigned int i; /* index of 1D array */
unsigned int j;
// draw critical point
ix = (int)((Zx-ZxMin)/PixelWidth);
iy = (int)((ZyMax-Zy)/PixelHeight); // reverse y axis
i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
data[i]=255-data[i];
// iterate
for (j = 1; j <= IterMax; j++) //larg number of iteration s
{ Zx2 = Zx*Zx;
Zy2 = Zy*Zy;
// bailout test
if (Zx2 + Zy2 > ER2) return iExterior; // if escaping stop iteration
// if not escaping iterate
// Z(n+1) = Zn * Zn + C
Zy = 2*Zx*Zy + Cy;
Zx = Zx2 - Zy2 + Cx;
//compute integer coordinate
ix = (int)((Zx-ZxMin)/PixelWidth);
iy = (int)((ZyMax-Zy)/PixelHeight); // reverse y axis
i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
data[i]=255-data[i]; // mark the trap
}
return 0;
}
// save data array to pgm file
/*
Bug :
"... when trying to compile your C code on my Mac off of
the wiki-books page, I consistently find that I need to increase the
string size for the file name. You always seem to allocate 10
characters for the file name and then use an sprintf to print a file
name to a string where the name is based on a float that you've
computed. On my Mac, that float always has more than 10 digits so,
when I try to run the program, I get an error. Increasing the buffer
size solved the problem." Mark McClure
*/
int SaveArray2PGMFile( unsigned char data[], double t)
{
FILE * fp;
const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ; it is 8 bit color file */
char name [10]; /* name of file */
sprintf(name,"%f", t); /* */
char *filename =strcat(name,".pgm");
char *comment="# ";/* comment should start with # */
/* save image to the pgm file */
fp= fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode */
fprintf(fp,"P5\n %s\n %u %u\n %u\n",comment,iWidth,iHeight,MaxColorComponentValue); /*write header to the file*/
fwrite(data,iSize,1,fp); /*write image data bytes to the file in one step */
printf("File %s saved. \n", filename);
fclose(fp);
return 0;
}
int info()
{
// diplay info messages
printf("InternalAngle = %f \n", InternalAngle);
printf("Cx = %f \n", Cx);
printf("Cy = %f \n", Cy);
//
printf("alfax = %f \n", creal(alfa));
printf("alfay = %f \n", cimag(alfa));
printf("iHeight = %d \n", iHeight);
printf("PixelWidth = %f \n", PixelWidth);
printf("distorsion of image = %f \n", ratio);
printf("iterMax = %lu \n", iterMax);
return 0;
}
/* ----------------------------------------- main -------------------------------------------------------------*/
int main()
{
setup();
// here are procedures for creating image file
//FillArray( data ); // no symmetry
FillArraySymmetric(data);
SaveArray2PGMFile(data , d2+0.000); // save edge array (only boundaries) to pgm file
AddBoundaries(data);
// CheckOrientation( data );
SaveArray2PGMFile(edge ,d2+0.001); // save edge array (only boundaries) to pgm file
CopyBoundaries();
SaveArray2PGMFile(data , d2+ 0.002); // save edge array (only boundaries) to pgm file
MarkTrap();
DrawCriticalOrbit(1000000);
SaveArray2PGMFile(data , d2+ 0.003); // save edge array (only boundaries) to pgm file
//
free(data);
free(edge);
//
info();
return 0;
}
Image Magic src code
convert 0.603000.pgm -resize 1000x1000 a.png
פריטים שמוצגים בקובץ הזה
מוצג
ערך כלשהו ללא פריט ויקינתונים
14 באפריל 2013
היסטוריית הקובץ
ניתן ללחוץ על תאריך/שעה כדי לראות את הקובץ כפי שנראה באותו זמן.
תאריך/שעה | תמונה ממוזערת | ממדים | משתמש | הערה | |
---|---|---|---|---|---|
נוכחית | 18:17, 14 באפריל 2013 | 1,000 × 1,000 (100 ק"ב) | Soul windsurfer | User created page with UploadWizard |
שימוש בקובץ
אין בוויקיפדיה דפים המשתמשים בקובץ זה.
שימוש גלובלי בקובץ
אתרי הוויקי השונים הבאים משתמשים בקובץ זה:
- שימוש באתר en.wikibooks.org
- שימוש באתר pl.wikipedia.org
מטא־נתונים
קובץ זה מכיל מידע נוסף, שכנראה הגיע ממצלמה דיגיטלית או מסורק שבהם הקובץ נוצר או עבר דיגיטציה.
אם הקובץ שונה ממצבו הראשוני, כמה מהנתונים להלן עלולים שלא לשקף באופן מלא את הקובץ הנוכחי.
הערה בקובץ PNG |
---|