/***********************************************
* Newtz6C.c
*
*  Generates a RAW format file of an image 
*  of a Julia set attractor of Newton's method (grayscale)
*
*  C source file by Alberto Strumia
*
************************************************/

#include<stdio.h>	

/* definition of constants */
#define Radius 0.7
#define Cx 0.0
#define Cy 0.0
#define Side 0.8
#define M 700
#define Num 128

/* main program */
main()
{

	int p, q, n, w;	
	double x, y, xx, yy, Incx, Incy;
	
	FILE *fp;
	
	fp = fopen("Newt.raw","w");

	for (p = 1; p <= M; p++)
	{
		Incy = - Side + 2*Side/M*p;
		
			printf("%i %%\n", p*100/M);	

		for (q = 1; q <= M; q++)
		{
			Incx = - Side + 2*Side/M*q;
			
		x = Incx;
		y = Incy;
		w = 0;
		
			for ( n = 1; n <= Num; ++n)
			{
			xx = 5*x/6.0 - x*(x*x*x*x - 10*x*x*y*y +5*y*y*y*y)/(x*x+y*y)/(x*x+y*y)/(x*x+y*y)/(x*x+y*y)/(x*x+y*y)/6.0;

			yy = 5*y/6.0 + y*(5*x*x*x*x - 10*x*x*y*y + y*y*y*y)/(x*x+y*y)/(x*x+y*y)/(x*x+y*y)/(x*x+y*y)/(x*x+y*y)/6.0;
				
				x = xx;
				y = yy;
				
				
				if ( (x-Cx)*(x-Cx) + (y-Cy)*(y-Cy) < Radius )	
					{
						w = n;
						n = Num;
					}
					
			}
		fprintf(fp, "%c", 2*w );
		}
	}
	fclose(fp);	
}

/* end of main program */