/* 
 c program:
  1. draws Mandelbrot set for Fm(z)=z^4+m*z;
  using  escape time ( boolean and levele sets )

Adam Majewski
fraktal.republka.pl
 
 -------------------------------         
 2. technic of creating ppm file is  based on the code of Claudio Rocchini
 http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
 create 24 bit color graphic file ,  portable pixmap file = PPM 
 see http://en.wikipedia.org/wiki/Portable_pixmap
 to see the file use external application ( graphic viewer)
 ---------------------------------

It is a console C program ( one file)
It can be compiled under :
*windows ( gcc thru Dev-C++ )
*linux and mac using gcc : 

 gcc m.c -lm -Wall
it creates a.out file. Then run it :
 ./a.out

It creates ppm file in program directory. Use file viewer to see it.







(%i1) z:zx+zy*%i;
(%o1) %i*zy+zx
(%i2) m:mx+my*%i;
(%o2) %i*my+mx
(%i3) z1:z^4+m*z;
(%o3) (%i*zy+zx)^4+(%i*my+mx)*(%i*zy+zx)
(%i4) realpart(z1);
(%o4) zy^4-6*zx^2*zy^2-my*zy+zx^4+mx*zx
(%i5) imagpart(z1);
(%o5) -4*zx*zy^3+4*zx^3*zy+mx*zy+my*zx
diff(z^4+m*z,z,1);


(%i1) diff(z^4+m*z,z,1);
(%o1) 4*z^3+m

s:to_poly_solve([4*z^3+m], [z]);
%union([z=-(%i*my+mx)^(1/3)/4^(1/3)],[z=-((sqrt(3)*%i-1)*(%i*my+mx)^(1/3))/(2*4^(1/3))],[z=((sqrt(3)*%i+1)*(%i*my+mx)^(1/3))/(2*4^(1/3))])

(%i8) b:first(s);
(%o8) [z=-(%i*my+mx)^(1/3)/4^(1/3)]

(%i18) b:rhs(b[1]);
(%o18) -(%i*my+mx)^(1/3)/4^(1/3)

(%i19) realpart(b);
(%o19) -((my^2+mx^2)^(1/6)*cos(atan2(my,mx)/3))/4^(1/3)
(%i20) imagpart(b);
(%o20) -((my^2+mx^2)^(1/6)*sin(atan2(my,mx)/3))/4^(1/3)



" z0 shall be a critical point, where the derivative is 0.
The derivative is  m + n*z^{n-1} so  z  is any {n-1}-th root of
-m/n . " Wolf Jung 

"
 I have made the changes
discussed below and it works.

1) 
In the iteration,  you have computed the new  Zy  and used it
to compute the new  Zx , where you should use the old  Zy .
Change this to
temp = ...
Zx = ...
Zy = temp

The iteration will be much faster,  if you compute  z^4  by
squaring  Z^2 :

for (Iteration=0; Iteration<IterationMax && Zx2+Zy2<ER2; Iteration++)
{  Zx2 -= Zy2; Zy2 = 2.0*Zx*Zy; //z^2
     double temp = 2.0*Zx2*Zy2 +Mx*Zy +My*Zx ;
     Zx = Zx2*Zx2 - Zy2*Zy2 - My*Zy + Mx*Zx;
     Zy = temp;
     Zx2=Zx*Zx;
     Zy2=Zy*Zy;
 };


2)
For the critical point it should be  atan2(-my, -mx)  and  Zy = +(b* ...
You have done complex conjugation instead of turning by 180°,
which -z means.

A minimal speed gain might be obtained by multiplying with the
cubic root of 1/4 instead of dividing by the cubic root of 4.

And I would write 1.0/6.0 instead of 0.16... 

Finally,  the cubic root of 4 could be in b:

a=atan2(-My,-Mx)/3.0; // atan2(-my,-mx)
b= pow(0.0625*(My*My + Mx*Mx) , 1.0/6.0);
Zx= b*cos(a);  Zy= b*sin(a);

   Best regards,

   Wolf ""






  */
 #include <stdio.h>
 #include <math.h>
 int main()
 {
          /* screen ( integer) coordinate */
        int iX,iY;
        const int iXmax = 5000; 
        const int iYmax = 5000;
        /* world ( double) coordinate = parameter plane*/
        double Mx,My;
        const double MxMin=-2.4;
        const double MxMax=2.4;
        const double MyMin=-2.4;
        const double MyMax=2.4;
        /* */
        double PixelWidth=(MxMax-MxMin)/iXmax;
        double PixelHeight=(MyMax-MyMin)/iYmax;
        /* color component ( R or G or B) is coded from 0 to 255 */
        /* it is 24 bit color RGB file */
        const int MaxColorComponentValue=255; 
        FILE * fp;
        char *filename="lsm50000.ppm";
        char *comment="# ";/* comment should start with # */
        static unsigned char color[3];
        /* Z=Zx+Zy*i  ;   Z0 = 0 */
        double Zx, Zy;
        double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
        /*  */
        int Iteration;
        const int IterationMax=10000;
        /* bail-out value , radius of circle ;  */
        const double EscapeRadius=3;
        double ER2=EscapeRadius*EscapeRadius;
        double a, b, temp;
        unsigned char ucTemp;
        
        /*create new file,give it a name and open it in binary mode  */
        fp= fopen(filename,"wb"); /* b -  binary mode */
        /*write ASCII header to the file*/
        fprintf(fp,"P6\n %s\n %d\n %d\n %d\n",comment,iXmax,iYmax,MaxColorComponentValue);
        
        /* compute and write image data bytes to the file*/
        for(iY=0;iY<iYmax;iY++)
        {
             My=MyMin + iY*PixelHeight;
             if (fabs(My)< PixelHeight/2) My=0.0; /* Main antenna */
             printf(" %d \n", iY);
             for(iX=0;iX<iXmax;iX++)
             {         
                        Mx=MxMin + iX*PixelWidth;
                        /* initial value of orbit = critical point Z:  (-m/4)^(1/3)=0  */
                        a=atan2(-My,-Mx)/3.0; // atan2(-my,-mx)
                        b= pow(0.0625*(My*My + Mx*Mx) , 1.0/6.0);
                        Zx=  b*cos(a);
                        Zy=  b*sin(a);
			//printf(" %f ; %f \n", Zx,Zy);
                        //printf("Zx = %f ; Zy = %f \n", Zx, Zy);
                        Zx2=Zx*Zx;
                        Zy2=Zy*Zy;
                        /* */
                        for (Iteration=0;Iteration<IterationMax && ((Zx2+Zy2)<ER2);Iteration++)
                        {
                            temp =-4.0*Zx*Zy*Zy2  +4.0*Zx2*Zx*Zy +Mx*Zy +My*Zx ; // -4*zx*zy^3 +4*zx^3*zy +mx*zy +my*zx
                            Zx=Zy2*Zy2 -6.0*Zx2*Zy2 +Zx2*Zx2 -My*Zy + Mx*Zx; // zy^4 -6*zx^2*zy^2 -my*zy +zx^4+ mx*zx
                            Zy=temp;
 			    //
			    Zx2=Zx*Zx;
                            Zy2=Zy*Zy;
                        };
                        /* compute  pixel color (24 bit = 3 bajts) */
                        if (Iteration<IterationMax)
                        { /* exterior of Mandelbrot set  */
                          ucTemp = 255 - 100*((unsigned char)(255*Iteration/IterationMax)) ;
                           color[0]=ucTemp; /* Red*/
                             color[1]=ucTemp;  /* Green */ 
                             color[2]=ucTemp;/* Blue */                          
                        }
                  
                        else /*  interior of Mandelbrot set  */
                           {
                             color[0]=0;
                           color[1]=0;
                           color[2]=0;     
                           };
                         /*write color to the file*/
                        fwrite(color,1,3,fp);
                }
        }
        fclose(fp);
        return 0;
 }