/* pointtovu: connect-up XYZ points into a vu file by distance.	 * Input:  lines with X,Y,Z in first three fields, or possibly	 * with an initial "s" or "l" field (for compatibility with .vu	 * files.) Lines beginning with a # or lines without three numeric	 * values are ignored (not even counted in line numbering)	 *  Mike Pique	 *	 * options:	 * -r radius     uses radius (float) as connection distance limit	 * Typical usage to make a bondlist might be: -r 1.85 	 *	 */#define DEFAULTDIST 1.8 /* connect points closer than this */#define MAXNPOINTS 14000#include <stdio.h>#include <math.h>int nflag = 1;int dflag = 0;int kflag = 0;int hflag = 0;double rad = DEFAULTDIST;#define HDIST 1.15char *pgmname;double atof();main(argc,argv)int argc;char **argv;{FILE *f = NULL; /* input file *//* option cracking: see getopt(3) */extern char *optarg;extern int optind;int c;	pgmname = argv[0];	rad = DEFAULTDIST;	while(( c = getopt(argc, argv, "r:ndkh")) != -1) switch(c) {	case 'r': rad = atof(optarg);		break;	case 'n': nflag++; break;	case 'd': dflag++; break;	case 'k': kflag++; break;	case 'h': hflag++; break;	case '?': usage(); break;	default: usage(); exit(-1); break;	}/*while for getopt*/	if(rad<0 || rad>1E9) {		fprintf(stderr,"%s: bad radius %g\n",		 pgmname, rad);		usage();		exit(-3);	}/*if radius no good*/	for(; optind < argc; optind++) {		f = fopen(argv[optind], "r");		if(f==NULL) {			fprintf(stderr,"%s: can't open %s.\n", pgmname, argv[optind]);			usage();			exit(-1);		}/*can't read file*/		do_dist(f, argv[optind]);		fclose(f);	}/*for each file argument*/	if(f == NULL) do_dist(stdin, "stdin");	exit(0);}/*main*/do_dist(f, filename)FILE * f;char *filename;{typedef struct {	float x, y, z;	char *key;	} point;static point pos[MAXNPOINTS];register point *point_ptr;float x, y, z;float dx, dy, dz;register float dist = rad;char linebuf[1000];/* char *malloc(); */register int p;int npoints, i;char *w;	npoints = 0;	while(fgets(linebuf,sizeof(linebuf),f) != NULL) {/*for each line of data in the input file do this*/		int nmatches = 0;		char keybuffer[1024];		float d; /* distance */		if(linebuf[0] == '#') continue;		if(linebuf[0] == 's' || linebuf[0]== 'l') linebuf[0]=' ';		if(3 > sscanf(linebuf, "%f %f %f %[^\n]", &x, &y, &z, keybuffer)) continue;		/* enter this point into list */		pos[npoints].x = x;		pos[npoints].y = y;		pos[npoints].z = z;		pos[npoints].key = (char *)malloc((size_t)(1+strlen(keybuffer)));	/* DEC Alpha */		/*pos[npoints].key = (char *)malloc((size_t)(1+strlen(keybuffer)));*/		if(pos[npoints].key != NULL) strcpy(pos[npoints].key, keybuffer);			#ifdef TRIPLET#define min(a,b) ((a)<(b)?(a):(b))#define max(a,b) ((a)>(b)?(a):(b))#define dist_pos(a,b) hypot(hypot(pos[a].x-pos[b].x, pos[a].y-pos[b].y),pos[a].z-pos[b].z)		{register int j,k;		for(j=0; j<npoints;j++) for(k=0;k<j;k++)  {			float dij, djk, dki;			float d1, d2, d3;			dij = dist_pos(npoints,j);			djk = dist_pos(j,k);			dki = dist_pos(npoints,k);			/* sort increasing - is there an easier way? mp */			d1 = min(dij,min(djk,dki)); /* min */			d2 = min(max(dij,djk), max(djk,dki)); /* middle */			d3 = max(dij,max(djk,dki)); /* max */			if(dflag)printf("%.4f %.4f %.4f", d1, d2, d3);#define out1(i) printf(" %g %g %g", pos[i].x, pos[i].y, pos[i].z)			out1(k);			out1(j);			out1(npoints);			if(kflag) printf(" %s %s %s",			  pos[k].key, pos[j].key, pos[npoints].key);			printf("\n");			}/*for j=0 to npoints*/		npoints++;/*now do whole loop again for one more point*/		continue;		}/*register j,k*/#endif /* TRIPLET */		/* search for points within 'd' of this one */		point_ptr = pos;		for( p= 0; p<npoints; p++ ) {			dx = x- point_ptr->x;			dy = y- point_ptr->y; 		w=pos[npoints].key+16;			if( (dx<dist|| dx> -dist) && (dy<dist|| dy> -dist) 			 && (d=hypot( z- point_ptr->z, hypot(dx,dy)))<dist &&			 (strncmp("H",point_ptr->key+16,1)||strncmp("H",pos[npoints].key+16,1))) {				/* accept */				if(nflag) printf("%d %d",p+1,npoints+1); /* 1-origin */				else				out(point_ptr->x,point_ptr->y, point_ptr->z,				 x,y,z);				if(kflag) printf(" %s %s", point_ptr->key, pos[npoints].key);				if(dflag) printf(" %g", d);				printf("\n");				nmatches++;				if(hflag&& d <= HDIST) break;				}			point_ptr++;/*printf("\n");*/			}		/* if no matches, emit the point as a point */		/* ?? what is right? mp  if(nmatches==0) out(x,y,z,x,y,z); */		npoints++;	}/*while fgets*/	if(kflag) for(p=0;p<npoints;p++) free(pos[p].key);}/*do_dist*/out(x1,y1,z1, x2, y2, z2)float x1,y1,z1,x2,y2,z2;{#define map(a) (a)	printf("%g %g %g %g %g %g",	map(x1), map(y1), map(z1),	map(x2), map(y2), map(z2));}/*out*/usage(){	fprintf(stderr,"%s: usage:  %s [-r radius (%g)] [-k] [-n] [-h]\n",	 pgmname, pgmname, DEFAULTDIST);}/*usage*/