#include <stdio.h>
#include <math.h>

int width=1024;
int height=768;
int aasamples=7;

const int ballcount = 44;
const int balls[][3] =
{
/*  0 */ { -15,   3,   9 },
/*  1 */ { -15,   0,   9 },
/*  2 */ { -15,  -3,   9 },
/*  3 */ { -15,  -6,   9 },
/*  4 */ {  -8,  -6,   9 },
/*  5 */ {  -8,   3,   9 },
/*  6 */ { -10,   0,   9 },
/*  7 */ {  -6,   0,   9 },
/*  8 */ { -10,  -3,   9 },
/*  9 */ {  -6,  -3,   9 },
/* 10 */ {  -1,   0,   9 },
/* 11 */ {  -1,  -3,   9 },
/* 12 */ {   1,  -6,   9 },
/* 13 */ {   1,   3,   9 },
/* 14 */ {   6,   0,   9 },
/* 15 */ {   6,  -3,   9 },
/* 16 */ {   8,  -6,   9 },
/* 17 */ {   8,   3,   9 },
/* 18 */ {  13,   0,   9 },
/* 19 */ {  13,  -3,   9 },
/* 20 */ {  15,  -6,   9 },
/* 21 */ {  15,   3,   9 },
/* 22 */ { -11, -11,   0 },
/* 23 */ { -11, -13,   0 },
/* 24 */ { -11, -15,   0 },
/* 25 */ { -11, -17,   0 },
/* 26 */ {  -9, -11,   0 },
/* 27 */ {  -9, -15,   0 },
/* 28 */ {  -7, -13,   0 },
/* 29 */ {  -7, -17,   0 },
/* 30 */ {   0, -11,   0 },
/* 31 */ {   0, -15,   0 },
/* 32 */ {  -2, -13,   0 },
/* 33 */ {  -2, -15,   0 },
/* 34 */ {  -2, -17,   0 },
/* 35 */ {   2, -13,   0 },
/* 36 */ {   2, -15,   0 },
/* 37 */ {   2, -17,   0 },
/* 38 */ {   7, -11,   0 },
/* 39 */ {   7, -13,   0 },
/* 40 */ {   9, -15,   0 },
/* 41 */ {   9, -17,   0 },
/* 42 */ {  11, -11,   0 },
/* 43 */ {  11, -13,   0 },
/* 44 */ { -64, -64,   0 },
};

float N=36;
float E;
float S;
float C;

void findball(int b)
{
	E = balls[b][0];
	C = balls[b][1];
	S = balls[b][2];
}

float bdist;

void q(int c, float x, float y, float z, float k, float l, float m, float a, float b)
{
	findball(c);
	x-=E;
	y-=S;
	z-=C;
	b =  x*x +  y*y +  z*z - 4.0;
	a = -x*k + -y*l + -z*m      ;
	b = a*a - b;
	if (b >= 0.0)
	{
		b = sqrt(b);
		if (a>b)
			bdist = a + -b;
		else
			bdist = a + b;
	} else
		bdist = -1.0;
}

float lastball;
float W;
void trace(float x, float y, float z, float k, float l, float m, int a)
{
	int ball;
	
	lastball=-1;
	for (ball=0; ball < ballcount; ball++)
	{
		q(ball,x,y,z,k,l,m,0,0);
		if (bdist>0 && ball!=a && (bdist<W || lastball<0))
		{
			W=bdist;
			lastball=ball;
		}
	}
}

float R;
float G;
float B;
float u;
float v;
float w;

void renderpxl(int iters, float x, float y, float z, float h, float i, float j, int a)
{
	float val;
	int ball;
	
	trace(x,y,z,h,i,j,a);
	if (iters < 9 && lastball>=0)
	{
		double distance;
		
		ball=lastball;
		findball(lastball);
		
		x+=h*W;
		y+=i*W;
		z+=j*W;
		
		u=x-E;
		v=y-S;
		w=z-C;
		
		val=(-2.0*u-2.0*v+w)/3.0;
		val/=2.0;
		val*=val;
		val*=200.0;
		
		distance = sqrt(u*u+v*v+w*w);
		if (distance!=0)
		{
			u=-u/distance;
			v=-v/distance;
			w=-w/distance;
		}
		
		E=(h*u+i*v+j*w);
		h-=u*E*2.0;
		i-=v*E*2.0;
		j-=w*E*2.0;
		renderpxl(iters+1,x,y,z,h,i,j,lastball);
		
		R/=2;
		G/=2;
		B/=2;
		
		     if (ball<22)  { 	R += val;	G += val;	B += val;	}
		else if (ball<30)  {	R += val;					}
		else if (ball<38)  {			G += val;			}
		else if (ball<44)  {					B += val;	}
		else if (ball==44) {			G += val;	B += val;	}
		else		   {	R += val;	G += val;			}
	} else if (!iters) {
		z+=2.0;
		if (z>0.0)
			j = z/8.0;
		else
			j = z/20.0;
		if (j>0)
		{
			val=j*j;
			R=255.0-250.0*val;
			G=255.0-150.0*val;
			B=255.0-100.0*val;
		} else {
			val=j*j;
			if (val < (1.0/5.0))
			{
				R=255.0-210.0*val;
				G=255.0-435.0*val;
				B=255.0-720.0*val;
			} else {
				val-=1.0/5.0;
				R=213.0-110.0*val;
				G=168.0-113.0*val;
				B=111.0-85.0*val;
			}
		}
	}
	
	if (R < 0) R = 0;
	if (R > 255) R = 255;
	
	if (G < 0) G = 0;
	if (G > 255) G = 255;
	
	if (B < 0) B = 0;
	if (B > 255) B = 255;
}

float J=0;
float K=-10;
float L=-7;


void render(int x, int y)
{
	int xaa, yaa;
	int raccum, gaccum, baccum;
	
	raccum = gaccum = baccum = 0;
	for (yaa = 0; yaa < aasamples; yaa++)
		for (xaa = 0; xaa < aasamples; xaa++)
		{
			renderpxl(
				0,
				40.0*(float)(aasamples*x+xaa)/(float)(aasamples*width )		- 20.0,
												- 10.0,
				-30.0*(float)(aasamples*y+yaa)/(float)(aasamples*height)	+ 8.0,
				0.0,
				1.0,
				0.0,
				-1);
			raccum+=R;
			gaccum+=G;
			baccum+=B;
		}
	raccum /= aasamples;	gaccum /= aasamples;	baccum /= aasamples;
	raccum /= aasamples;	gaccum /= aasamples;	baccum /= aasamples;
	printf("%c%c%c",raccum,gaccum,baccum);
}

int main()
{
	int y,x;
	printf("P6\n%i %i\n255\n", width, height);
	for (y = 0; y < height; y++)
	{
		for (x=0; x < width; x++)
		{
			if ((x % 17) == 0)	// make it look semi-random
			{
				fprintf(stderr, "Rendering: row %4d, col %4d   \r", y, x);
				fflush(stderr);
			}
			render(x,y);
		}
	}
	fprintf(stderr, "Done                                                \n");
	return 0;
}