/* Modified 13 feb 1999
 * Modified 23 feb 1999 (neatened)
 * Modified 6 oct 2000 (added scaling)
 * How to run me:
Compile using
  'cc mandel.c'
then type:
  'a.out'
a.out will send the image via STDOUT, so you'll need to redirect its
output (e.g. '>' ) to a file.

*/

#include <stdio.h>
#include <math.h>
/*#include <rasterfile.h>*/
#define RAS_MAGIC 0x59a66a95L
#define RT_OLD 0
#define RMT_NONE 0


#define ULONGMAX 0x7fffffffL
/*#define PI 3.14159265*/



float rnd();
void writeraster();
void prl();

long seed1=422316127,seed2=25393472;
int nr_of_iter=255,image_x_size=10240/2,image_y_size=7680/2;
int scale_x=5, scale_y=5;
float *r;



main(argv,argc)
int argv;
char *argc[];
{
int a,b,c,d;
double x,y,nx,ny, aa,bb;
double x_phase,y_phase, x_mult,y_mult;
double top,left,xsize,ysize;
/* Not working. These are from an older program.
  if (argv>1) sscanf(argc[1],"%d",&seed1);
  if (argv>2) sscanf(argc[2],"%d",&seed2);
  if (argv>3) sscanf(argc[3],"%d",&nr_of_iter);
  if (argv>4) sscanf(argc[4],"%d",&image_x_size);
*/
	/* add guard lines around image (these are not saved) */
	image_x_size += scale_x*2;
	image_y_size += scale_y*2;
	r=(float *)calloc((image_y_size/scale_y)*(image_x_size/scale_x),
			sizeof(float) );

 	/*** mandel1:
	top=-0.25;
	left=-0.75;
	size=0.05;
	***/  
    /*** mandel2: 
	left=-0.5760719498822918;
	top=-0.629;
	xsize=-0.0033;
    ***/
	/*** mandel3:
	left=0.095;
	top=0.622;
	xsize=0.12-left;
	***/
	/*** mandel4:
	left=-1.269;
	top=0.0398392165219373;
	xsize=0.01;
	***/
	/*** mandel5:
	***/
	left = 0.3095;
	top=-0.0296;
	xsize=0.005;
	/*** mandel6:
	left=-0.7497;
	top=-0.131;
	xsize=0.012;
	***/
	/*** mandel7:
	left=-0.1759075362495006;
	top=0.6492722801549806;
	xsize=0.0012;
	***/
	/*** mandel8:
	left=-0.8637;
	top=0.2640841426373636;
	xsize=0.002;
	***/
	/*** mandel9
	left=-1.4027558440885592;
	top=0.0;
	xsize=0.0025;
	***/
    /*** mandel10
	left=-0.42;
	top=0.594;
	xsize=0.016;
	***/

	ysize=xsize * (double)image_y_size / (double)image_x_size;

	for (b=0; b<image_y_size - scale_y; b++) {
		for (a=0; a<image_x_size - scale_x; a++) {
			aa=left + xsize  / (double)image_x_size*(double)a;
			bb=top + ysize  / (double)image_y_size*(double)b;
			x=aa;
			y=bb;
			for (c=nr_of_iter;c>1;c--) {
				nx=x*x -y*y + aa;
				ny=2*x*y + bb;
				if (nx*nx + ny*ny > 4) break;
				x=nx;
				y=ny;
			}
		/*if (c=0) c=255;*/
			/*if (a<scale_x || a>=image_x_size - 2*scale_x) {
				x_mult=0.25;
			} else {
				x_mult=0.25;
			}
			if (b<scale_y || b>=image_y_size - 2*scale_y) {
				y_mult=0.25;
			} else {
				y_mult=0.25;
			}*/

			x_mult=0.25;
			y_mult=0.25;
	
			x_phase = x_mult*((double)(a % scale_x))/scale_x;
			y_phase = y_mult*((double)(b % scale_y))/scale_y;
			r[ a         /scale_x+(image_x_size/scale_x)*( b         /scale_y)]
			  += (x_mult-x_phase + y_mult-y_phase)*c;
			r[(a+scale_x)/scale_x+(image_x_size/scale_x)*( b         /scale_y)]
			  += (       x_phase + y_mult-y_phase)*c;
			r[ a         /scale_x+(image_x_size/scale_x)*((b+scale_y)/scale_y)]
			  += (x_mult-x_phase +        y_phase)*c;
			r[(a+scale_x)/scale_x+(image_x_size/scale_x)*((b+scale_y)/scale_y)]
			  += (       x_phase +        y_phase)*c;
		}
	}

	writeraster();
}





void writeraster() {
long i,j;
/*SUNraster (for more information see 'man rasterfile') */
/*Header*/
  prl(RAS_MAGIC);
  prl(image_x_size/scale_x - 2);
  prl(image_y_size/scale_y - 2);
  prl(8);
  prl((image_x_size/scale_x - 2)*(image_y_size/scale_y - 2));
  prl(RT_OLD);
  prl(RMT_NONE);
  prl(0); 
/*cmap:standard*/

/*image:*/

  for (i=1; i<image_y_size/scale_y - 1; i++)
    for (j=1; j<image_x_size/scale_x - 1; j++)
      printf("%c",(int)r[(image_x_size/scale_x)*i + j]/(scale_x*scale_y));

}

	void
prl (long l) {
	char c;
	printf("%c%c%c%c",(l>>24)&255, (l>>16) & 255, (l>>8) & 255, l& 255);
}


float rnd(mul)
  float mul;
{
  seed2=(seed2>>2)^(seed1<<2)+0xd576493b;
  seed1-=seed2+0x93a46b53;
  return ( ((float)(seed2 & ULONGMAX)) / ULONGMAX * mul );
}