/*******************************************************************

rotate.c -- Eric R. Weeks, started 4/4/97

v01:  4-04-97: taken from shell.c, rotate.gwk

This program is public domain, although please leave my name,
email, & web info attached.

weeks@physics.emory.edu
http://www.physics.emory.edu/~weeks/software/rotate.html

 *******************************************************************/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define PI 3.14159265358979323846264338328
#define EE 2.71828182845904523536
#define MAX 222222

char *prgname;

main(argc,argv)
int argc;
char **argv;
{
	int c,numpts,t,r;
	double x,y,z,xx,yy,zz;
	double max[3],min[3],center[3];
	FILE *fp;
	extern int optind;
	extern char *optarg;
	float data[MAX][3];
	float theta,phi,eta;
	float ct,st,cp,sp,ce,se;
	int centerofmass,usebox,spacing;

	centerofmass = 1;			/* rotate around center of mass  */
	usebox = 1;				/* do use a box                  */
	spacing = 40;				/* number of points per box edge */
	theta = 0.0;
	phi = 0.0;
	eta = 0.0;
	while ((c = getopt(argc, argv, "hp:t:CBs:r:")) != EOF)
	switch (c)  {
	case 'h': fprintf(stderr,"Usage: %s [options] [<] [file]\n\n",argv[0]);
fprintf(stderr," -h   : this help message\n");
fprintf(stderr,"\nview from this angle (theta measured from z-axis):\n");
fprintf(stderr," -t # : read in theta value (degrees) [%.1f]\n",theta);
fprintf(stderr," -p # : read in phi value (degrees) [%.1f]\n",phi);
fprintf(stderr," -r # : third rotation angle (degrees) [%.1f]\n",eta);
fprintf(stderr," -C   : do NOT rotate around center of mass\n");
fprintf(stderr," -B   : do NOT draw a box around data\n");
fprintf(stderr," -s # : change spacing of box points [%d]\n",spacing);
		   exit(1);
		   break;
	case 't': theta = atof(optarg);
			theta = theta * PI/180.0;
			break;
	case 'p': phi = atof(optarg);
			phi = phi * PI/180.0;
			break;
	case 'C': centerofmass = 0;
			break;
	case 'B': usebox = 0;
			break;
	case 's': spacing = atoi(optarg);
			break;
	case 'r': eta = atof(optarg);
			eta = eta * PI / 180.0;
			break;
        }

	argc -= (optind-1) ; argv += (optind-1) ;
	fp = (argc > 1) ? fopen(*++argv, "r") : stdin;

	numpts = 0;
	for (r=0;r<3;r++)  {
		center[r] = 0.0;
		max[r] = -9e99;
		min[r] = 9e99;
	}
	while ((fscanf(fp,"%lf %lf %lf",&x,&y,&z)) == 3)  {
		data[numpts][0] = x;
		data[numpts][1] = y;
		data[numpts][2] = z;
		for (r=0;r<3;r++)  {
			center[r]+=data[numpts][r];
			if (data[numpts][r]>max[r])  max[r] = data[numpts][r];
			if (data[numpts][r]<min[r])  min[r] = data[numpts][r];
		}
		++numpts;
	}
	if (numpts<1)  {
		fprintf(stderr,"No data read in, exiting...\n");
		exit(1);
	}
	if (centerofmass==1)
		for (r=0;r<3;r++)  center[r] /= ((float)numpts);
	else
		for (r=0;r<3;r++)  center[r] = 0.0;

	if (usebox==1)  {
		for (t=0;t<spacing;t++)  {
			x = ((float)t)/((float)spacing);
			xx = min[0]; yy = min[1]; zz = min[2]*x + max[2]*(1.0-x);
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;
			xx = max[0]; yy = min[1]; zz = min[2]*x + max[2]*(1.0-x);
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;
			xx = max[0]; yy = max[1]; zz = min[2]*x + max[2]*(1.0-x);
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;
			xx = min[0]; yy = max[1]; zz = min[2]*x + max[2]*(1.0-x);
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;

			xx = min[0]; yy = min[1]*x + max[1]*(1.0-x); zz = min[2];
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;
			xx = max[0]; yy = min[1]*x + max[1]*(1.0-x); zz = min[2];
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;
			xx = max[0]; yy = min[1]*x + max[1]*(1.0-x); zz = max[2];
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;
			xx = min[0]; yy = min[1]*x + max[1]*(1.0-x); zz = max[2];
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;

			xx = min[0]*x + max[0]*(1.0-x); yy = min[1]; zz = min[2];
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;
			xx = min[0]*x + max[0]*(1.0-x); yy = min[1]; zz = max[2];
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;
			xx = min[0]*x + max[0]*(1.0-x); yy = max[1]; zz = max[2];
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;
			xx = min[0]*x + max[0]*(1.0-x); yy = max[1]; zz = min[2];
			data[numpts][0]=xx;  data[numpts][1]=yy;  data[numpts][2]=zz;
			numpts++;
		}
	}

	ct = cos(theta);
	st = sin(theta);
	cp = cos(phi);
	sp = sin(phi);
	ce = cos(eta);
	se = sin(eta);

	for (t=0;t<numpts;t++)  {
		/* printf("%f %f %f\n",data[t][0],data[t][1],data[t][2]); */
		/* printf("%f %f %f\n", center[0], center[1], center[2]); */
		/* printf("%f %f %f %f %f %f\n",theta,ct,st,phi,cp,sp); */
		xx = data[t][0]-center[0];
		yy = data[t][1]-center[1];
		zz = data[t][2]-center[2];
		x = ct*xx - st*zz;  y = yy;            z = st*xx + ct*zz;
		xx = x;             yy = cp*y - sp*z;  zz = sp*y + cp*z;
		x = ce*xx - se*yy;  y = se*xx + ce*yy; z = zz;
		xx = x + center[0];
		yy = y + center[1];
		zz = z + center[2];
		printf("%f %f %f\n",xx,yy,zz);
	}
}