/* Sphere Inversion sample generator (VRML 1.0)
   (C)2012 Claudio Rocchini */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

const double PI = 3.1415926535897932384626433832795;

class point3 {
public:
	double x,y,z;
	point3() {}
	point3( double nx, double ny, double nz ): x(nx),y(ny),z(nz) {}
};
double norm( const point3 & p ) {
	return sqrt(p.x*p.x + p.y*p.y + p.z*p.z);
}
point3 operator+ ( const point3 & a, const point3 & b ) {
	return point3( a.x+b.x, a.y+b.y, a.z+b.z );
}
point3 operator- ( const point3 & a, const point3 & b ) {
	return point3( a.x-b.x, a.y-b.y, a.z-b.z );
}
point3 operator* ( const point3 & a, double s ) {
	return point3( a.x*s, a.y*s, a.z*s );
}

point3 inversion( const point3 & CS, double RS, const point3 & p ) {
	point3 d = p-CS;
	double dp = norm(d);
	double k = (RS*RS)/(dp*dp);
	return CS + d*k;
}

int main() {
	const double GS = 0.1;				// global scale
	point3 CS(0,0,0); double RS = 100;	// sphere center and radius
	point3 CC(0,0,0); point3 AC(0,1,0);	// cylinder "center" and axis
	point3 A1(1,0,0); point3 A2(0,0,1);	// cylinder hortobase
	double RC = 50;	double HC = 1600;	// cyinder radius and height
	
	FILE * fo = fopen("sphereinv.wrl","w");
	fprintf(fo,"#VRML V1.0 ascii\n");
	fprintf(fo,
		"Separator {\n"
		"  Material {\n"
		"    ambientColor   0.2 0.2 0.2\n diffuseColor   0.4 0.9 0.2\n"
		"    specularColor  0.3 0.3 0.3\n emissiveColor  0 0 0\n"
		"    shininess      0.4\n transparency   0.5\n"
		"  }\n"
		"  Translation { translation %g %g %g }\n"
		"  Sphere { radius %g }\n"
		"}\n"
		,CS.x*GS,CS.y*GS,CS.z*GS,RS*GS
	);
	fprintf(fo,"Coordinate3 {\n""point [\n");
	
	const int N = 32;
	const int M = 256;
	int i,j;
	for(i=0;i<N;++i)
	for(j=0;j<M;++j)
	{
		double p1 = 2*PI*i/N;		// 1st parameter
		double c1 = cos(p1)*RC;
		double s1 = sin(p1)*RC;
		double p2 = -HC/2 + HC*j/M;	// 2nd parameter
		point3 p = CC + A1*c1 + A2*s1 + AC*p2;
		fprintf(fo,"%g %g %g,\n",p.x*GS,p.y*GS,p.z*GS);
		point3 ip = inversion(CS,RS,p);
		fprintf(fo,"%g %g %g,\n",ip.x*GS,ip.y*GS,ip.z*GS);
	}
	fprintf(fo,"]\n}\n");
	fprintf(fo,"Separator {\n");
	fprintf(fo,
		"  Material {\n"
		"    ambientColor   0 0 0\n diffuseColor   0 0 0\n"
		"    specularColor  0 0 0\n emissiveColor  0 0 0.8\n"
		"    shininess      0.0\n transparency   0.0\n"
		"  }\n"
		);
	fprintf(fo,"IndexedLineSet {\ncoordIndex [\n");
	for(i=0;i<N;++i) for(j=0;j<M-1;++j) fprintf(fo,"%d,%d,-1,\n",(j+i*M)*2,((j+1)+i*M)*2);
	for(j=0;j<M;++j) for(i=0;i<N;++i)   fprintf(fo,"%d,%d,-1,\n",(j+i*M)*2,(j+((i+1)%N)*M)*2);
	fprintf(fo,"]\n}\n}\nSeparator {\n");
	fprintf(fo,
		"  Material {\n"
		"    ambientColor   0 0 0\n diffuseColor   0 0 0\n"
		"    specularColor  0 0 0\n emissiveColor  0.8 0 0\n"
		"    shininess      0.0\n transparency   0.0\n"
		"  }\n"
		);
	fprintf(fo, "IndexedLineSet {\n coordIndex [\n");
	for(i=0;i<N;++i) for(j=0;j<M;++j) fprintf(fo,"%d,%d,-1,\n",(j+i*M)*2+1,((j+1)%M+i*M)*2+1);
	for(j=0;j<M;++j) for(i=0;i<N;++i) fprintf(fo,"%d,%d,-1,\n",(j+i*M)*2+1,(j+((i+1)%N)*M)*2+1);
	fprintf(fo, "]\n}\n}\n");
	fclose(fo);
	return 0;
}