/***********************************************
* Julia3Dsilv.c
*
*  Generates a RAW format file of a 3D image 
*  of a Julia set (grayscale)
*
*  C source file by Alberto Strumia
*
************************************************/


#include<stdio.h>

/* definition of constants */

#define Radius 10	
#define Side 1.2	
#define Cx 0.745429		
#define Cy 0.05	
#define scale 2		
#define M 700		
#define Num 127		
#define co 1.41	
#define si 1.41		

/* main program */
main()
{

	int p, q, n, wext, wint, r[M+Num/scale];	
	double x, y, xx, yy, Incx, Incy, md;	

	FILE *fp;
	fp = fopen("Julia.raw","w");	
	
	for (p = 1; p <= M; p++)
	{
		
			printf("%i %%\n", p*100/M);	
		
		for (q = 1; q <= M; q++)
		{
			Incx = -2*Side + 2*Side*(p*co+q*si/2)/M;
			Incy = 0.5*Side + 2*Side*(-p*si+q*co/2)/M;
			
		x = Incx;
		y = Incy;

			for ( n = 1; n <= Num; ++n)
			{
				xx = x*x - y*y - Cx;
				yy = 2*x*y - Cy;
				
				x = xx;
				y = yy;
				md = x*x + y*y;
				
				wint = md*Num*10;

			if ( wint < Num )
			     {
				wext = wint;
			      }
			else
			 {
			     wint = 150;
			 }

				r[q] = wint;

				if ( md > Radius )
				  {
					wext = n;
					n = Num;
				  }

			}
		
		for ( n = 1; n <= wext; n++)
			{
				if ( q > n/scale )
				   {
					r[q + n/scale] = wext;
				   }
			}
		}
				
		for ( q = 1; q <= M; q++)
			{	
			
		fprintf(fp, "%c", 2*r[q] );		
			}
	}
	
	fclose(fp);			

/* end of main program */