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

#define MAXOBJ 10
#define EPSILON (0.005)
#define SIN60 (0.866025404)
#define RADIUS (0.5 - EPSILON)			// avoid contact between spheres
#define DELTA (1.0 / (2.0 * SIN60))
#define SUBRADIUS (DELTA - RADIUS - 2*EPSILON)

typedef struct spherestruct {
	double x;
	double y;
	double z;
	double radius;
} Sphere;

Sphere sphere[MAXOBJ];

/*
** Output a height field of hexagonally-packed hemispheres.
*/
int main(int argc, char** argv) {
	int i, j;
	int obj;
	int resolution = 8;
	int repfactor = 1;
	int objects = MAXOBJ;
	float vertscale = 1.0;
	double x, y;
	double dist;
	double height;
	double maxheight;

	/*
	** Evaluate command line arguments
	*/
	for (i=0; i < argc; ++i) {
		if ((argv[i][0] == '-') && (argv[i][1] == 'r')) {
			resolution = atoi(&argv[i][2]);
		}
		else if ((argv[i][0] == '-') && (argv[i][1] == 'f')) {
			repfactor = atoi(&argv[i][2]);
		}
		else if ((argv[i][0] == '-') && (argv[i][1] == 'v')) {
			sscanf(&argv[i][2], "%f", &vertscale);
		}
		else if ((argv[i][0] == '-') && (argv[i][1] == 'o')) {
			objects = atoi(&argv[i][2]);
		}
	}
	fprintf(stderr, "resolution = %d\n", resolution);
	fprintf(stderr, "objects = %d\n", objects);
	fprintf(stderr, "replication factor = %d\n", repfactor);
	fprintf(stderr, "vertical scale = %f\n", vertscale);

	/*
	** Initialize the sphere structures
	**
	** Pattern is hexagonal mesh of spheres of radius 1/2, tightly packed
	**   (these alone cover 90% of the surface)
	** Smaller spheres are located in the intersticial spaces
	**   (these cover an addition 5% of the surface)
	*/
	sphere[0].x = 0.5;
	sphere[0].y = 0.0;
	sphere[0].z = 0.0;
	sphere[0].radius = RADIUS;
	sphere[1].x = 0.0;
	sphere[1].y = SIN60;
	sphere[1].z = 0.0;
	sphere[1].radius = RADIUS;
	sphere[2].x = 1.0;
	sphere[2].y = SIN60;
	sphere[2].z = 0.0;
	sphere[2].radius = RADIUS;
	sphere[3].x = 0.5;
	sphere[3].y = 2 * SIN60;
	sphere[3].z = 0.0;
	sphere[3].radius = RADIUS;

	sphere[4].x = 0.0;
	sphere[4].y = SIN60 - DELTA;
	sphere[4].z = 0.0;
	sphere[4].radius = SUBRADIUS;
	sphere[5].x = 1.0;
	sphere[5].y = SIN60 - DELTA;
	sphere[5].z = 0.0;
	sphere[5].radius = SUBRADIUS;
	sphere[6].x = 0.5;
	sphere[6].y =  DELTA;
	sphere[6].z = 0.0;
	sphere[6].radius = SUBRADIUS;
	sphere[7].x = 0.5;
	sphere[7].y = 2 * SIN60 - DELTA;
	sphere[7].z = 0.0;
	sphere[7].radius = SUBRADIUS;
	sphere[8].x = 0.0;
	sphere[8].y = SIN60 + DELTA;
	sphere[8].z = 0.0;
	sphere[8].radius = SUBRADIUS;
	sphere[9].x = 1.0;
	sphere[9].y = SIN60 + DELTA;
	sphere[9].z = 0.0;
	sphere[9].radius = SUBRADIUS;

	/*
	** Step through the entire mesh
	*/
	for (j=0; j < (resolution * repfactor); ++j) {
		y = 2.0 * SIN60 * (double)(j % resolution) / resolution;
		for (i=0; i < (resolution * repfactor); ++i) {
			x = (double)(i % resolution) / resolution;
			maxheight = 0.0;
			for (obj=0; obj < objects; ++obj) {
				dist = sqrt((x - sphere[obj].x) *
					    (x - sphere[obj].x) +
					    (y - sphere[obj].y) *
					    (y - sphere[obj].y));
				if (dist <= sphere[obj].radius) {
					height = sphere[obj].radius *
						sin(acos(dist /
							sphere[obj].radius));
					if (maxheight < height)
						maxheight = height;
				}
			}
		printf("%8.6f ", maxheight * vertscale);
		}
	printf("\n");
	}
	return 0;
}