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

/* geodome.c 
 * 
 * (C) 1993 Andy Wardley (abw@oasis.icl.co.uk).  All Rights Reserved.
 *
 * Creates triangle descriptions for geodomes for POVray v2 ray-tracer.
 * Optional command line parameter specifes number of sub-divisions to
 * make (default is 3).  Requires math library to compile on some machines.
 * 
 * Source and/or executables may be freely distributed.  Any part of this
 * code may be re-hashed or used as long as credit is given.  No guarantee
 * is provided - take it as you find it.  POVray files created may also
 * be distributed, etc.  Any comments, POV files or general ray-tracing
 * chit-chat - email me at above address.
 *
 * Enjoy.
 * A
 */ 


struct vertex
{
	double x, y, z;
};

struct triangle
{
	struct vertex a, b, c;	
};


struct vertex mid_point(struct vertex *v1, struct vertex *v2)
{
	struct vertex midv;
	double ratio;

	/* find the mid-point by averaging the two ends */
	midv.x = (v1->x + v2->x) / 2;
	midv.y = (v1->y + v2->y) / 2;
	midv.z = (v1->z + v2->z) / 2;

	/* now "stretch" the point out so that the distance from origin is 1 */
	ratio = sqrt(pow(midv.x, 2) + pow(midv.y, 2) + pow(midv.z, 2));
	midv.x /= ratio;
	midv.y /= ratio;
	midv.z /= ratio;

	return midv;
}


void divide_triangle(struct triangle *t, int depth)
{
	struct triangle new_t[4];
	struct vertex m1, m2, m3;
	int nt;

	/* assign the outside corners of the triangle to the sub-triangles */
	new_t[0].a = t->a;
	new_t[1].b = t->b;
	new_t[2].c = t->c;

	/* calculate and assign the mid-points */
	new_t[0].b = new_t[1].a = new_t[3].c = mid_point(&t->a, &t->b);
	new_t[0].c = new_t[2].a = new_t[3].b = mid_point(&t->a, &t->c);
	new_t[1].c = new_t[2].b = new_t[3].a = mid_point(&t->b, &t->c);

	if (--depth == 0)	/* depth reached - print */
	{
		for(nt = 0; nt < 4; nt++)
		{
			printf("\ttriangle {\n");
			printf("\t\t<%.10f, %.10f, %.10f>\n",
				new_t[nt].a.x, new_t[nt].a.y, new_t[nt].a.z);
			printf("\t\t<%.10f, %.10f, %.10f>\n",
				new_t[nt].b.x, new_t[nt].b.y, new_t[nt].b.z);
			printf("\t\t<%.10f, %.10f, %.10f>\n",
				new_t[nt].c.x, new_t[nt].c.y, new_t[nt].c.z);
			printf("\t}\n");
		}
	}
	else			/* recurse and sub-divide */
		for(nt = 0; nt < 4; nt++)
			divide_triangle(&new_t[nt], depth);
}


void main(int argc, char **argv)
{
	int depth;
	struct triangle t;

	if (argc < 2 || (depth = atoi(argv[1])) == 0)
		depth = 3;

	/* define the corners of a triangle <0, 1, 0> <0, 0, 1> <1, 0, 0> */
	t.a.x = t.b.x = 0;
	t.b.y = t.c.y = 0;
	t.a.z = t.c.z = 0;
	t.a.y = t.b.z = t.c.x = 1;

	/* print object header */
	printf("/* triangle data created by the geodome utility */\n");
	printf("/* writen by Andy Wardley (abw@oasis.icl.co.uk) */\n\n");
	printf("#declare geo_part =\nunion {\n");

	/* calculate the triangles */
	divide_triangle(&t, depth);

	/* print object trailer */
	printf("}\n\n#declare geodome =\nobject {\n\tunion {\n");
	printf("\t\tobject { geo_part }\n");
	printf("\t\tobject { geo_part rotate <0, 90, 0> }\n");
	printf("\t\tobject { geo_part rotate <0, 180, 0> }\n");
	printf("\t\tobject { geo_part rotate <0, 270, 0> }\n");
	printf("\t\tobject { geo_part rotate <180, 0, 0> }\n");
	printf("\t\tobject { geo_part rotate <180, 90, 0> }\n");
	printf("\t\tobject { geo_part rotate <180, 180, 0> }\n");
	printf("\t\tobject { geo_part rotate <180, 270, 0> }\n");
	printf("\t}\n}\n\n");
}