Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- <?php
- $text = "== {{int:filedesc}} ==\r\n{{Information\r\n|Description ={{en|1=Parabolic Julia set for internal angle 1 over 5 }}\r\n|Source =own work with help of : see comments in code \r\n|Author =[[User:Adam majewski|Adam majewski]]\r\n|Date =2012-10-07\r\n|Permission =\r\n|other_versions =\r\n}}\r\n\r\n==C src code ==\r\n<source lang=c>\r\n/*\r\n\r\nAdam Majewski\r\nfraktal.republika.pl\r\n\r\nc console progam using \r\n* symmetry\r\n* openMP\r\n\r\nIt uses modified DEM method to draw parabolic julia sets\r\n\r\n\r\n\r\ngcc t.c -lm -Wall -fopenmp -march=native \r\ntime ./a.out\r\n\r\n\r\nFile s1000000f1.5.pgm saved. \r\nCx = 0.356763 \r\nCy = 0.328582 \r\nalfax = 0.154508 \r\nalfay = 0.475528 \r\ndistorsion of image = 1.000000 \r\n\r\nreal\t123m19.106s\r\n----------------------------\r\nFile s10000000f1.5.pgm saved. \r\nCx = 0.356763 \r\nCy = 0.328582 \r\nalfax = 0.154508 \r\nalfay = 0.475528 \r\ndistorsion of image = 1.000000 \r\n\r\nreal\t1023m48.596s = 17 fours\r\n-------------------------------------\r\nso s100000000f1.5.pgm shoud take 17 days !!!!!! \r\n\r\n\r\n\r\n\r\n\r\n\r\n*/\r\n\r\n\r\n\r\n#include <stdio.h>\r\n#include <stdlib.h> // malloc\r\n#include <string.h> // strcat\r\n#include <math.h> // M_PI; needs -lm also \r\n#include <complex.h>\r\n#include <omp.h> // OpenMP; needs also -fopenmp\r\n\r\n\r\n/* --------------------------------- global variables and consts ------------------------------------------------------------ */\r\n\r\n\r\n// virtual 2D array and integer ( screen) coordinate\r\n// Indexes of array starts from 0 not 1 \r\nunsigned int ix, iy; // var\r\nunsigned int ixMin = 0; // Indexes of array starts from 0 not 1\r\nunsigned int ixMax ; //\r\nunsigned int iWidth ; // horizontal dimension of array\r\nunsigned int ixAxisOfSymmetry ; // \r\nunsigned int iyMin = 0; // Indexes of array starts from 0 not 1\r\nunsigned int iyMax ; //\r\nunsigned int iyAxisOfSymmetry ; // \r\nunsigned int iyAbove ; // var, measured from 1 to (iyAboveAxisLength -1)\r\nunsigned int iyAboveMin = 1 ; //\r\nunsigned int iyAboveMax ; //\r\nunsigned int iyAboveAxisLength ; //\r\nunsigned int iyBelowAxisLength ; //\r\nunsigned int iHeight = 1000; // odd number !!!!!! = (iyMax -iyMin + 1) = iyAboveAxisLength + iyBelowAxisLength +1\r\n// The size of array has to be a positive constant integer \r\nunsigned int iSize ; // = iWidth*iHeight; \r\n\r\n\r\n// memmory 1D array \r\nunsigned char *data;\r\n// unsigned int i; // var = index of 1D array\r\nunsigned int iMin = 0; // Indexes of array starts from 0 not 1\r\nunsigned int iMax ; // = i2Dsize-1 = \r\n// The size of array has to be a positive constant integer \r\n// unsigned int i1Dsize ; // = i2Dsize = (iMax -iMin + 1) = ; 1D array with the same size as 2D array\r\n\r\n\r\n/* world ( double) coordinate = dynamic plane */\r\n const double ZxMin=-1.5;\r\n const double ZxMax=1.5;\r\n const double ZyMin=-1.5;\r\n const double ZyMax=1.5;\r\n double PixelWidth; // =(ZxMax-ZxMin)/iXmax;\r\n double PixelHeight; // =(ZyMax-ZyMin)/iYmax;\r\n double distanceMax;\r\n double ratio ;\r\n double lambda=1.5;\r\n\r\n\r\n\r\n// complex numbers of parametr plane \r\ndouble Cx; // c =Cx +Cy * i\r\ndouble Cy;\r\ndouble complex c; // \r\n\r\ndouble complex alfa; // alfa fixed point alfa=f(alfa)\r\n\r\nunsigned long int iterMax = 10000000; //iHeight*100;\r\ndouble ER = 2.0; // Escape Radius for bailout test \r\ndouble ER2;\r\n\r\n\r\n/* ------------------------------------------ functions -------------------------------------------------------------*/\r\n\r\n\r\n\r\n/* find c in component of Mandelbrot set \r\n \r\n uses code by Wolf Jung from program Mandel\r\n see function mndlbrot::bifurcate from mandelbrot.cpp\r\n http://www.mndynamics.com/indexp.html\r\n\r\n */\r\ndouble complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int period)\r\n{\r\n //0 <= InternalRay<= 1\r\n //0 <= InternalAngleInTurns <=1\r\n double t = InternalAngleInTurns *2*M_PI; // from turns to radians\r\n double R2 = InternalRadius * InternalRadius;\r\n //double Cx, Cy; /* C = Cx+Cy*i */\r\n switch ( period ) // of component \r\n {\r\n case 1: // main cardioid\r\n Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4; \r\n Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4; \r\n break;\r\n case 2: // only one component \r\n Cx = InternalRadius * 0.25*cos(t) - 1.0;\r\n Cy = InternalRadius * 0.25*sin(t); \r\n break;\r\n // for each period there are 2^(period-1) roots. \r\n default: // higher periods : to do\r\n Cx = 0.0;\r\n Cy = 0.0; \r\n break; }\r\n\r\n return Cx + Cy*I;\r\n}\r\n\r\n\r\n/*\r\n\r\nhttp://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings\r\nz^2 + c = z\r\nz^2 - z + c = 0\r\nax^2 +bx + c =0 // ge3neral for of quadratic equation\r\nso :\r\na=1\r\nb =-1\r\nc = c\r\nso :\r\n\r\nThe discriminant is the d=b^2- 4ac \r\n\r\nd=1-4c = dx+dy*i\r\nr(d)=sqrt(dx^2 + dy^2)\r\nsqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx +- sy*i\r\n\r\nx1=(1+sqrt(d))/2 = beta = (1+sx+sy*i)/2\r\n\r\nx2=(1-sqrt(d))/2 = alfa = (1-sx -sy*i)/2\r\n\r\nalfa : attracting when c is in main cardioid of Mandelbrot set, then it is in interior of Filled-in Julia set, \r\nit means belongs to Fatou set ( strictly to basin of attraction of finite fixed point )\r\n\r\n*/\r\n// uses global variables : \r\n// ax, ay (output = alfa(c)) \r\ndouble complex GiveAlfaFixedPoint(double complex c)\r\n{\r\n double dx, dy; //The discriminant is the d=b^2- 4ac = dx+dy*i\r\n double r; // r(d)=sqrt(dx^2 + dy^2)\r\n double sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i\r\n double ax, ay;\r\n \r\n // d=1-4c = dx+dy*i\r\n dx = 1 - 4*creal(c);\r\n dy = -4 * cimag(c);\r\n // r(d)=sqrt(dx^2 + dy^2)\r\n r = sqrt(dx*dx + dy*dy);\r\n //sqrt(d) = s =sx +sy*i\r\n sx = sqrt((r+dx)/2);\r\n sy = sqrt((r-dx)/2);\r\n // alfa = ax +ay*i = (1-sqrt(d))/2 = (1-sx + sy*i)/2\r\n ax = 0.5 - sx/2.0;\r\n ay = sy/2.0;\r\n \r\n\r\nreturn ax+ay*I;\r\n}\r\n\r\n\r\ndouble GiveDistanceBetween(double complex z1, double complex z2 )\r\n{double dx,dy;\r\n \r\n dx = creal(z1) - creal(z2);\r\n dy = cimag(z1) - cimag(z2);\r\n return sqrt(dx*dx+dy*dy);\r\n \r\n} \r\n\r\n\r\n\r\nint setup()\r\n{\r\n\r\n unsigned int period = 5; // period of secondary component joined by root point\r\n unsigned int denominator;\r\n double InternalAngle;\r\n \r\n\r\n denominator = period;\r\n InternalAngle = 1.0/((double) denominator);\r\n\r\n c = GiveC(InternalAngle, 1.0, 1) ;\r\n Cx=creal(c);\r\n Cy=cimag(c);\r\n alfa = GiveAlfaFixedPoint(c);\r\n\r\n /* 2D array ranges */\r\n if (!(iHeight % 2)) iHeight+=1; // it sholud be even number (variable % 2) or (variable & 1)\r\n iWidth = iHeight;\r\n iSize = iWidth*iHeight; // size = number of points in array \r\n // iy\r\n iyMax = iHeight - 1 ; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].\r\n iyAboveAxisLength = (iHeight -1)/2;\r\n iyAboveMax = iyAboveAxisLength ; \r\n iyBelowAxisLength = iyAboveAxisLength; // the same \r\n iyAxisOfSymmetry = iyMin + iyBelowAxisLength ; \r\n // ix\r\n \r\n ixMax = iWidth - 1;\r\n\r\n /* 1D array ranges */\r\n // i1Dsize = i2Dsize; // 1D array with the same size as 2D array\r\n iMax = iSize-1; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].\r\n\r\n\r\n /* Pixel sizes */\r\n PixelWidth = (ZxMax-ZxMin)/ixMax; // ixMax = (iWidth-1) step between pixels in world coordinate \r\n PixelHeight = (ZyMax-ZyMin)/iyMax;\r\n ratio = ((ZxMax-ZxMin)/(ZyMax-ZyMin))/((float)iWidth/(float)iHeight); // it should be 1.000 ...\r\n distanceMax = PixelWidth; \r\n //\r\n \r\n ER2 = ER * ER;\r\n\r\n /* create dynamic 1D arrays for colors ( shades of gray ) */\r\n \r\n data = malloc( iSize * sizeof(unsigned char) );\r\n if (data == NULL )\r\n {\r\n fprintf(stderr,\" Could not allocate memory\\n\");\r\n getchar();\r\n return 1;\r\n }\r\n else fprintf(stderr,\" memory is OK \\n\");\r\n\r\n \r\n \r\n return 0;\r\n\r\n}\r\n\r\n\r\n\r\n// from screen to world coordinate ; linear mapping\r\n// uses global cons\r\ndouble GiveZx(unsigned int ix)\r\n{ return (ZxMin + ix*PixelWidth );}\r\n\r\n// uses globaal cons\r\ndouble GiveZy(unsigned int iy)\r\n{ return (ZyMax - iy*PixelHeight);} // reverse y axis\r\n\r\n\r\n// forward iteration of complex quadratic polynomial\r\n// fc(z) = z*z +c\r\n// z0 = critical point \r\n// uses global var \r\nlong int GiveLastIteration(double Zx, double Zy )\r\n{ \r\n long int iter;\r\n // z= Zx + Zy*i\r\n double Zx2, Zy2; /* Zx2=Zx*Zx; Zy2=Zy*Zy */\r\n \r\n /* initial value of orbit = critical point Z= 0 */\r\n \r\n Zx2=Zx*Zx;\r\n Zy2=Zy*Zy;\r\n // for iter from 0 to iterMax because of ++ after last loop\r\n for (iter=0; iter<iterMax && ((Zx2+Zy2)<ER2); ++iter)\r\n {\r\n Zy=2*Zx*Zy + Cy;\r\n Zx=Zx2-Zy2 + Cx;\r\n Zx2=Zx*Zx;\r\n Zy2=Zy*Zy;\r\n };\r\n \r\n return iter;\r\n}\r\n\r\n\r\n\r\n\r\n \r\n\r\n\r\n\r\n/*\r\n estimates distance from point c to nearest point in Julia set \r\n for Fc(z)= z*z + c\r\n z(n+1) = Fc(zn) \r\n this function is based on function mndlbrot::dist from mndlbrot.cpp\r\n from program mandel by Wolf Jung (GNU GPL )\r\n http://www.mndynamics.com/indexp.html \r\n\r\nHyunsuk Kim : \r\nFor Julia sets, z is the variable and c is a constant. Therefore df[n+1](z)/dz = 2*f[n]*f'[n] -- you don't add 1.\r\n-----------------\r\n\r\n\r\n\"The following method was taught to me by Christian Henriksen.\r\nIt is a distance estimator method (different from Milnor's).\r\nWhen iterating forward, keep also tract of the value of the derivative.\r\nMore precisely when computing z_n from z_0, also compute dz_n/dz_0.\r\nby the chain rule, it is equal to the product of the derivatives of f taken at z_0, z_1, .. , up to z_n-1\r\nTherefore you can compute it recursively alongside and like the value of z_n.\r\nThen choose a point in the Julia set, here I chose the critical point.\r\nThen if the ball :\r\n- of center z_n \r\n- and radius lambda * size of a pixel * derivative \r\ncontains this point there is a good chance \r\nthat the pixel of center z_0 \r\ncontains an iterated preimage of the critical point, thus I color z_0 in black.\r\n\r\nHere lambda is a thickness factor that you must choose. \r\nnot too small otherwise you do not see enough points, but not too big because otherwise you get artifacts.\r\n\r\nBest Regards,\r\nArnaud.\"\r\n\r\n\r\n */\r\n int jdist(double Zx, double Zy, double Cx, double Cy , int iter_max)\r\n { \r\n int i;\r\n double x = Zx, /* Z = x+y*i */\r\n y = Zy, \r\n /* Zp = xp+yp*1 = 1 */\r\n xp = 1, \r\n yp = 0, \r\n /* temporary */\r\n nz, \r\n nzp,\r\n /* a = abs(z) */\r\n //a,\r\n radius,\r\n distance; \r\n\r\n\r\n \r\n for (i = 1; i <= iter_max; i++)\r\n { /* first derivative zp = 2*z*zp = xp + yp*i; */\r\n \r\n nz = 2*(x*xp - y*yp) ; \r\n yp = 2*(x*yp + y*xp); \r\n xp = nz;\r\n /* z = z*z + c = x+y*i */\r\n nz = x*x - y*y + Cx; \r\n y = 2*x*y + Cy; \r\n x = nz; \r\n //\r\n nz = (x*x + y*y); // abs2(z)\r\n nzp = (xp*xp + yp*yp);\r\n if ( nz > ER2) return 0; // if escapes //nzp > 1e60 ||\r\n if (nzp>1e60) return 0; // ? interior ?\r\n distance=GiveDistanceBetween(alfa,x+y*I);\r\n nzp = sqrt(xp*xp + yp*yp);\r\n radius=lambda*PixelWidth*nzp;\r\n \r\n if (distance<radius)return 1; // \r\n }\r\n \r\n return 0; // \r\n }\r\n\r\n\r\nunsigned char GiveColor(unsigned int ix, unsigned int iy)\r\n{ \r\n double Zx, Zy; // Z= Zx+ZY*i;\r\n unsigned char color; // gray from 0 to 255 \r\n// unsigned char ColorList[]={255,230,180}; /* shades of gray used in imaage\r\n //color=0;\r\n // from screen to world coordinate \r\n Zx = GiveZx(ix);\r\n Zy = GiveZy(iy);\r\n //if ( GiveLastIteration(Zx, Zy ) == iterMax)\r\n // color = 0; // interior\r\n // else { \r\n\r\n// only dist\r\n if (jdist(Zx, Zy, Cx, Cy , iterMax))\r\n color = 0; // black = boundary = Julia set\r\n else color = 255; // Fatou set\r\n //}\r\n return color;\r\n}\r\n\r\n\r\n/* ----------- array functions -------------- */\r\n\r\n\r\n/* gives position of 2D point (iX,iY) in 1D array ; uses also global variable iWidth */\r\nunsigned int Give_i(unsigned int ix, unsigned int iy)\r\n{ return ix + iy*iWidth; }\r\n// ix = i % iWidth;\r\n// iy = (i- ix) / iWidth;\r\n// i = Give_i(ix, iy);\r\n\r\n\r\n\r\n\r\n// plots raster point (ix,iy) \r\nint PlotPoint(unsigned int ix, unsigned int iy, unsigned char iColor)\r\n{\r\n unsigned i; /* index of 1D array */\r\n i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */\r\n data[i] = iColor;\r\n\r\nreturn 0;\r\n}\r\n\r\n\r\n// fill array \r\n// uses global var : ...\r\n// scanning complex plane \r\nint FillArray(unsigned char data[] )\r\n{\r\n unsigned int ix, iy; // pixel coordinate \r\n\r\n\r\n// for all pixels of image \r\nfor(iy = iyMin; iy<=iyMax; ++iy) \r\n { printf(\" %d z %d\\n\", iy, iyMax); //info \r\n for(ix= ixMin; ix<=ixMax; ++ix) PlotPoint(ix, iy, GiveColor(ix, iy) ); // \r\n } \r\n \r\n return 0;\r\n}\r\n\r\n\r\n// fill array using symmetry of image \r\n// uses global var : ...\r\nint FillArraySymmetric(unsigned char data[] )\r\n{\r\n \r\n unsigned char Color; // gray from 0 to 255 \r\n\r\nprintf(\"axis of symmetry \\n\"); \r\niy = iyAxisOfSymmetry; \r\n#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)\r\nfor(ix=ixMin;ix<=ixMax;++ix) {//printf(\" %d from %d\\n\", ix, ixMax); //info \r\n PlotPoint(ix, iy, GiveColor(ix, iy));\r\n}\r\n\r\n\r\n/*\r\nThe use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.\r\nThe use of ‘private(variable, variable2)’ specifies that these variables should have a seperate instance in each thread.\r\n*/\r\n\r\n#pragma omp parallel for schedule(dynamic) private(iyAbove,ix,iy,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)\r\n\r\n// above and below axis \r\nfor(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) \r\n {printf(\" %d from %d\\n\", iyAbove, iyAboveMax); //info \r\n for(ix=ixMin; ix<=ixMax; ++ix) \r\n\r\n { // above axis compute color and save it to the array\r\n iy = iyAxisOfSymmetry + iyAbove;\r\n Color = GiveColor(ix, iy);\r\n PlotPoint(ix, iy, Color ); \r\n // below the axis only copy Color the same as above without computing it \r\n PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color ); \r\n } \r\n} \r\n return 0;\r\n}\r\n\r\n\r\n\r\n\r\n// Check Orientation of image : first quadrant in upper right position\r\n// uses global var : ...\r\nint CheckOrientation(unsigned char data[] )\r\n{\r\n unsigned int ix, iy; // pixel coordinate \r\n double Zx, Zy; // Z= Zx+ZY*i;\r\n unsigned i; /* index of 1D array */\r\n for(iy=iyMin;iy<=iyMax;++iy) \r\n {\r\n Zy = GiveZy(iy);\r\n for(ix=ixMin;ix<=ixMax;++ix) \r\n {\r\n\r\n // from screen to world coordinate \r\n Zx = GiveZx(ix);\r\n i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */\r\n if (Zx>0 && Zy>0) data[i]=255-data[i]; // check the orientation of Z-plane by marking first quadrant */\r\n\r\n }\r\n }\r\n \r\n return 0;\r\n}\r\n\r\n\r\n\r\n// save data array to pgm file \r\nint SaveArray2PGMFile( unsigned char data[])\r\n{\r\n \r\n FILE * fp;\r\n const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ; it is 8 bit color file */\r\n char name [10]; /* name of file */\r\n sprintf(name,\"s%luf%3.1f\", iterMax, lambda); /* */\r\n char *filename =strcat(name,\".pgm\");\r\n char *comment=\"# \";/* comment should start with # */\r\n\r\n /* save image to the pgm file */ \r\n fp= fopen(filename,\"wb\"); /*create new file,give it a name and open it in binary mode */\r\n fprintf(fp,\"P5\\n %s\\n %u %u\\n %u\\n\",comment,iWidth,iHeight,MaxColorComponentValue); /*write header to the file*/\r\n fwrite(data,iSize,1,fp); /*write image data bytes to the file in one step */\r\n printf(\"File %s saved. \\n\", filename);\r\n fclose(fp);\r\n\r\n return 0;\r\n}\r\n\r\n\r\nint info()\r\n{\r\n // diplay info messages\r\n printf(\"Cx = %f \\n\", Cx); \r\n printf(\"Cy = %f \\n\", Cy);\r\n // \r\n \r\n printf(\"alfax = %f \\n\", creal(alfa));\r\n printf(\"alfay = %f \\n\", cimag(alfa));\r\n printf(\"distorsion of image = %f \\n\", ratio);\r\nreturn 0;\r\n}\r\n\r\n\r\n/* ----------------------------------------- main -------------------------------------------------------------*/\r\nint main()\r\n{\r\n// here are procedures for creating image file\r\n\r\n setup();\r\n // compute colors of pixels = image\r\n //FillArray( data ); // no symmetry\r\n FillArraySymmetric(data); \r\n // CheckOrientation( data );\r\n SaveArray2PGMFile( data); // save array (image) to pgm file \r\n free(data);\r\n info();\r\n return 0;\r\n}\r\n</source>\r\n\r\n[[Category:Uploaded with UploadWizard]]\r\n[[Category:Parabolic Julia sets]]\r\n[[Category:Images with C source code]]\r\n[[Category:Black and white]]\r\n[[Category:Complex quadratic map]]\r\n\r\n== {{int:license-header}} ==\r\n{{self|cc-by-sa-3.0|GFDL}}.";
- echo "\r\n\r\n***Outputting the text in order to show that this is a valid string***\r\n";
- echo $text;
- echo "\r\n\r\n***Outputting a regular expression run on the text, which returns null***\r\n";
- var_dump(preg_replace("/([\s\S]+?\s*?)\s*\s*(\{\{)/u", "Hello world", $text));
- ?>
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement