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;
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);
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)
422 for (k=
dim-count; k--; )
424 }
else if (count >
dim) {
430 while ((s= strchr(s,
'P'))) {
438 if (issimplex+issimplex2) {
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) {
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);
562 norm= sqrt(1.0 +
x*
x +
y*
y);
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) {
598 anglediff= 2*
PI/numpoints;
603 offset= sqrt(radius * radius - (1-gap)*(1-gap)) - lensbase;
604 for (i=0; i < numpoints; i++, angle += anglediff) {
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) {
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) {
797 return (
int)(
a - 0.5);
799 if (
a + 0.5 > INT_MAX) {
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;