28 #pragma warning( disable : 4706) 29 #pragma warning( disable : 4996) 33 #define PI 3.1415926535897932384 39 void qh_out3n(
double a,
double b,
double c);
40 void qh_outcoord(
int iscdd,
double *coord,
int dim);
41 void qh_outcoincident(
int coincidentpoints,
double radius,
int iscdd,
double *coord,
int dim);
93 int coincidentcount=0, coincidenttotal=0, coincidentpoints=0;
94 int cubesize, diamondsize,
seed=0, count, apex;
95 int dim=3, numpoints=0, totpoints, addpoints=0;
96 int issphere=0, isaxis=0, iscdd=0, islens=0, isregular=0, iswidth=0, addcube=0;
97 int isgap=0, isspiral=0, NOcommand=0,
adddiamond=0;
98 int israndom=0, istime=0;
99 int isbox=0, issimplex=0, issimplex2=0, ismesh=0;
100 double width=0.0, gap=0.0, radius=0.0, coincidentradius=0.0;
101 double coord[
MAXdim], offset, meshm=3.0, meshn=4.0, meshr=5.0;
102 double *coordp, *simplex= NULL, *simplexp;
103 int nthroot, mult[
MAXdim];
104 double norm, factor, randr, rangap, lensangle=0, lensbase=1;
105 double anglediff, angle,
x,
y, cube=0.0, diamond=0.0;
108 char command[200], seedbuf[200];
109 char *s= command, *
t, *first_point= NULL;
114 qh_fprintf_rbox(rbox.
ferr, 6188,
"rbox error: rbox in use by another process. Please lock calls to rbox.\n");
121 exitcode= setjmp(rbox.
errexit);
131 strncat(command, rbox_command,
sizeof(command)-strlen(command)-1);
133 while (*s && !isspace(*s))
136 while (*s && isspace(*s))
211 if (*s && !isspace(*s)) {
212 qh_fprintf_rbox(rbox.
ferr, 7080,
"rbox error: arguments for 'Cn,r,m' are not 'int', 'float', and 'int'. Remaining string is '%s'\n", s);
215 if (coincidentpoints==0){
216 qh_fprintf_rbox(rbox.
ferr, 6268,
"rbox error: missing arguments for 'Cn,r,m' where n is the number of coincident points, r is the radius (default 0.0), and m is the number of points\n");
219 if (coincidentpoints<0 || coincidenttotal<0 || coincidentradius<0.0){
220 qh_fprintf_rbox(rbox.
ferr, 6269,
"rbox error: negative arguments for 'Cn,m,r' where n (%d) is the number of coincident points, m (%d) is the number of points, and r (%.2g) is the radius (default 0.0)\n", coincidentpoints, coincidenttotal, coincidentradius);
259 meshr= sqrt(meshn*meshn + meshm*meshm);
260 if (*s && !isspace(*s)) {
261 qh_fprintf_rbox(rbox.
ferr, 7069,
"rbox warning: assuming 'M3,4,5' since mesh args are not integers or reals\n");
262 meshn= 3.0, meshm=4.0, meshr=5.0;
272 while (*s && !isspace(*s))
287 qh_fprintf_rbox(rbox.
ferr, 7070,
"rbox error: unknown flag at %s.\nExecute 'rbox' without arguments for documentation.\n", s);
290 if (*s && !isspace(*s)) {
300 cubesize= (int)floor(ldexp(1.0,dim)+0.5);
317 qh_fprintf_rbox(rbox.
ferr, 6191,
"rbox error: lens radius %.2g should be greater than 1.0\n",
321 lensangle= asin(1.0/radius);
322 lensbase= radius * cos(lensangle);
328 else if (isregular + issimplex + islens + issphere + isaxis + isspiral + iswidth + ismesh) {
338 if ((issimplex + islens + isspiral + ismesh > 1)
339 || (issimplex + issphere + isspiral + ismesh > 1)) {
340 qh_fprintf_rbox(rbox.
ferr, 6193,
"rbox error: can only specify one of 'l', 's', 'x', 'Ln', or 'Mn,m,r' ('Ln s' is ok).\n");
343 if (coincidentpoints>0 && (numpoints == 0 || coincidenttotal > numpoints)) {
344 qh_fprintf_rbox(rbox.
ferr, 6270,
"rbox error: 'Cn,r,m' requested n coincident points for each of m points. Either there is no points or m (%d) is greater than the number of points (%d).\n", coincidenttotal, numpoints);
347 if (coincidenttotal == 0)
348 coincidenttotal= numpoints;
351 if (issimplex || ismesh)
352 totpoints= numpoints;
354 totpoints= numpoints+dim+1;
355 else if (isregular) {
356 totpoints= numpoints;
359 totpoints += numpoints - 2;
360 }
else if (dim == 3) {
362 totpoints += 2 * numpoints;
364 totpoints += 1 + numpoints;
369 totpoints= numpoints + isaxis;
370 totpoints += cubesize + diamondsize + addpoints;
371 totpoints += coincidentpoints*coincidenttotal;
375 for (s=command; *s; s++) {
376 if (issimplex2 && *s ==
'y')
382 }
else if (israndom) {
383 seed= (int)time(&timedata);
384 sprintf(seedbuf,
" t%d", seed);
385 strncat(command, seedbuf,
sizeof(command)-strlen(command)-1);
386 t= strstr(command,
" t ");
396 NOcommand ?
"" : command,
406 if ((s= first_point)) {
414 if (isspace(*s) || !*s)
417 qh_fprintf_rbox(rbox.
ferr, 6194,
"rbox error: missing comma after coordinate in %s\n\n", s);
422 for (k=dim-count; k--; )
424 }
else if (count > dim) {
425 qh_fprintf_rbox(rbox.
ferr, 6195,
"rbox error: %d coordinates instead of %d coordinates in %s\n\n",
430 while ((s= strchr(s,
'P'))) {
438 if (issimplex+issimplex2) {
439 if (!(simplex= (
double*)
qh_malloc( dim * (dim+1) *
sizeof(
double)))) {
445 for (i=0; i<dim; i++) {
446 for (k=0; k<dim; k++)
447 *(simplexp++)= i==k ? 1.0 : 0.0;
449 for (k=0; k<dim; k++)
452 for (i=0; i<dim+1; i++) {
453 for (k=0; k<dim; k++) {
455 *(simplexp++)= 2.0 * randr/randmax - 1.0;
461 for (i=0; i<dim+1; i++) {
464 for (k=0; k<dim; k++)
469 for (j=0; j<numpoints; j++) {
474 for (k=0; k<dim; k++)
477 for (i=0; i<dim+1; i++) {
479 factor= randr/randmax;
483 for (k=0; k<dim; k++) {
484 simplexp= simplex + i*dim + k;
485 coord[k] += factor * (*simplexp);
488 for (k=0; k<dim; k++)
489 coord[k] *= box/norm;
491 if(coincidentcount++ < coincidenttotal)
500 nthroot= (int)(pow((
double)numpoints, 1.0/dim) + 0.99999);
503 for (i=0; i < numpoints; i++) {
505 for (k=0; k < dim; k++) {
507 *(coordp++)= mult[0] * meshn + mult[1] * (-meshm);
509 *(coordp++)= mult[0] * meshm + mult[1] * meshn;
511 *(coordp++)= mult[k] * meshr;
514 if(coincidentcount++ < coincidenttotal)
516 for (k=0; k < dim; k++) {
517 if (++mult[k] < nthroot)
524 else if (isregular && !islens) {
525 if (dim != 2 && dim != 3) {
527 qh_fprintf_rbox(rbox.
ferr, 6197,
"rbox error: regular points can be used only in 2-d and 3-d\n\n");
530 if (!isaxis || radius == 0.0) {
545 anglediff= 2.0 *
PI/numpoints;
546 for (i=0; i < numpoints; i++) {
548 x= radius * cos(angle);
549 y= radius * sin(angle);
555 norm= sqrt(1.0 + x*x + y*y);
558 qh_out3n( box*x/norm, box*y/norm, box/norm);
562 norm= sqrt(1.0 + x*x + y*y);
565 qh_out3n( box*x/norm, box*y/norm, box/norm);
571 else if (isregular && islens && dim == 2) {
575 anglediff= 2 * lensangle/(numpoints - 1);
576 cos_0= cos(lensangle);
577 for (i=0; i < numpoints; i++, angle -= anglediff) {
578 x= radius * sin(angle);
579 y= radius * (cos(angle) - cos_0);
583 if (i != 0 && i != numpoints - 1) {
591 else if (isregular && islens && dim != 2) {
594 qh_fprintf_rbox(rbox.
ferr, 6198,
"rbox error: regular points can be used only in 2-d and 3-d\n\n");
598 anglediff= 2*
PI/numpoints;
603 offset= sqrt(radius * radius - (1-gap)*(1-gap)) - lensbase;
604 for (i=0; i < numpoints; i++, angle += anglediff) {
614 qh_out3n( box*x, box*y, box * offset);
617 qh_out3n( box*x, box*y, -box * offset);
626 for (j=0; j < gendim; j++)
635 for (i=0; i < numpoints; i++) {
637 for (j=0; j < gendim; j++) {
639 coord[j]= 2.0 * randr/randmax - 1.0;
640 norm += coord[j] * coord[j];
650 rangap= 1.0 - gap * randr/randmax;
651 factor= radius * rangap / norm;
652 for (j=0; j<gendim; j++)
653 coord[j]= factor * coord[j];
655 }
else if (islens && issphere) {
661 rangap= 1.0 - gap * randr/randmax;
662 factor= rangap / norm;
663 for (j=0; j<gendim; j++)
664 coord[j]= factor * coord[j];
666 }
else if (islens && !issphere) {
673 coord[j]= -1.0 - coord[j] * gap;
675 coord[j]= 1.0 - coord[j] * gap;
677 }
else if (isspiral) {
680 qh_fprintf_rbox(rbox.
ferr, 6199,
"rbox error: spiral distribution is available only in 3d\n\n");
683 coord[0]= cos(2*
PI*i/(numpoints - 1));
684 coord[1]= sin(2*
PI*i/(numpoints - 1));
685 coord[2]= 2.0*(double)i/(
double)(numpoints-1) - 1.0;
687 }
else if (issphere) {
691 factor *= 1.0 - width * randr/randmax;
693 for (j=0; j<dim; j++)
694 coord[j]= factor * coord[j];
697 if (isaxis && issphere) {
700 for (j=0; j<gendim; j++)
701 norm += coord[j] * coord[j];
703 for (j=0; j<dim; j++)
704 coord[j]= coord[j] / norm;
707 coord[dim-1] *= 1 - width * randr/randmax;
710 }
else if (isaxis && !issphere) {
712 coord[dim-1]= 2.0 * randr/randmax - 1.0;
715 coord[dim-1]= lensbase;
716 for (j=0, norm= 0; j<dim; j++)
717 norm += coord[j] * coord[j];
719 for (j=0; j<dim; j++)
720 coord[j]= coord[j] * radius/ norm;
721 coord[dim-1] -= lensbase;
724 coord[dim-1] *= 1 - width * randr/randmax;
727 coord[dim-1]= -coord[dim-1];
729 }
else if (iswidth && !issphere) {
732 coord[j]= -1.0 - coord[j] * width;
734 coord[j]= 1.0 - coord[j] * width;
737 for (k=0; k<dim; k++)
738 coord[k]= coord[k] * box;
742 if(coincidentcount++ < coincidenttotal)
749 for (j=0; j<cubesize; j++) {
752 for (k=dim-1; k>=0; k--) {
764 for (j=0; j<diamondsize; j++) {
767 for (k=dim-1; k>=0; k--) {
793 if (a - 0.5 < INT_MIN) {
794 qh_fprintf_rbox(rbox.
ferr, 6200,
"rbox input error: negative coordinate %2.2g is too large. Reduce 'Bn'\n", a);
797 return (
int)(a - 0.5);
799 if (a + 0.5 > INT_MAX) {
800 qh_fprintf_rbox(rbox.
ferr, 6201,
"rbox input error: coordinate %2.2g is too large. Reduce 'Bn'\n", a);
803 return (
int)(a + 0.5);
837 for (k=0; k < dim; k++)
842 void qh_outcoincident(
int coincidentpoints,
double radius,
int iscdd,
double *coord,
int dim) {
848 for (i= 0; i<coincidentpoints; i++) {
852 for (k=0; k < dim; k++) {
854 delta= 2.0 * randr/randmax - 1.0;
868 longjmp(rbox.
errexit, exitcode);
void * qh_malloc(size_t size)
void seed(unsigned int seed_value)
double qh_strtod(const char *s, char **endp)
void qh_fprintf_rbox(FILE *fp, int msgcode, const char *fmt,...)
void qh_outcoord(int iscdd, double *coord, int dim)
void qh_errexit_rbox(int exitcode)
void qh_out3n(double a, double b, double c)
int qh_strtol(const char *s, char **endp)
int qh_rboxpoints(FILE *fout, FILE *ferr, char *rbox_command)
void qh_out2n(double a, double b)
void qh_outcoincident(int coincidentpoints, double radius, int iscdd, double *coord, int dim)
void adddiamond(coordT *points, int numpoints, int numnew, int dim)
#define qh_RANDOMseed_(seed)