28 #pragma warning( disable : 4706)
29 #pragma warning( disable : 4996)
33 #define PI 3.1415926535897932384
41 void qh_outcoincident(
qhT *
qh,
int coincidentpoints,
double radius,
int iscdd,
double *coord,
int dim);
75 int coincidentcount=0, coincidenttotal=0, coincidentpoints=0;
76 int cubesize, diamondsize,
seed=0, count, apex;
77 int dim=3 , numpoints= 0, totpoints, addpoints=0;
78 int issphere=0, isaxis=0, iscdd= 0, islens= 0, isregular=0, iswidth=0, addcube=0;
79 int isgap=0, isspiral=0, NOcommand= 0,
adddiamond=0;
80 int israndom=0, istime=0;
81 int isbox=0, issimplex=0, issimplex2=0, ismesh=0;
82 double width=0.0, gap=0.0, radius=0.0, coincidentradius=0.0;
83 double coord[
MAXdim], offset, meshm=3.0, meshn=4.0, meshr=5.0;
84 double *coordp, *simplex= NULL, *simplexp;
86 double norm, factor, randr, rangap, lensangle= 0, lensbase= 1;
87 double anglediff, angle,
x,
y, cube= 0.0, diamond= 0.0;
90 char command[200], seedbuf[200];
91 char *s= command, *
t, *first_point= NULL;
95 exitcode= setjmp(
qh->rbox_errexit);
104 strncat(command, rbox_command,
sizeof(command)-strlen(command)-1);
106 while (*s && !isspace(*s))
109 while (*s && isspace(*s))
167 qh->rbox_isinteger= 1;
184 if (*s && !isspace(*s)) {
185 qh_fprintf_rbox(
qh,
qh->ferr, 7080,
"rbox error: arguments for 'Cn,r,m' are not 'int', 'float', and 'int'. Remaining string is '%s'\n", s);
188 if (coincidentpoints==0){
189 qh_fprintf_rbox(
qh,
qh->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");
192 if (coincidentpoints<0 || coincidenttotal<0 || coincidentradius<0.0){
193 qh_fprintf_rbox(
qh,
qh->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);
232 meshr= sqrt(meshn*meshn + meshm*meshm);
233 if (*s && !isspace(*s)) {
234 qh_fprintf_rbox(
qh,
qh->ferr, 7069,
"rbox warning: assuming 'M3,4,5' since mesh args are not integers or reals\n");
235 meshn= 3.0, meshm=4.0, meshr=5.0;
245 while (*s && !isspace(*s))
260 qh_fprintf_rbox(
qh,
qh->ferr, 7070,
"rbox error: unknown flag at %s.\nExecute 'rbox' without arguments for documentation.\n", s);
263 if (*s && !isspace(*s)) {
264 qh_fprintf_rbox(
qh,
qh->ferr, 7071,
"rbox error: missing space between flags at %s.\n", s);
270 if (
qh->rbox_isinteger && !isbox)
273 cubesize= (int)floor(ldexp(1.0,
dim)+0.5);
290 qh_fprintf_rbox(
qh,
qh->ferr, 6191,
"rbox error: lens radius %.2g should be greater than 1.0\n",
294 lensangle= asin(1.0/radius);
295 lensbase= radius * cos(lensangle);
301 else if (isregular + issimplex + islens + issphere + isaxis + isspiral + iswidth + ismesh) {
311 if ((issimplex + islens + isspiral + ismesh > 1)
312 || (issimplex + issphere + isspiral + ismesh > 1)) {
313 qh_fprintf_rbox(
qh,
qh->ferr, 6193,
"rbox error: can only specify one of 'l', 's', 'x', 'Ln', or 'Mn,m,r' ('Ln s' is ok).\n");
316 if (coincidentpoints>0 && (numpoints == 0 || coincidenttotal > numpoints)) {
317 qh_fprintf_rbox(
qh,
qh->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);
320 if (coincidenttotal == 0)
321 coincidenttotal= numpoints;
324 if (issimplex || ismesh)
325 totpoints= numpoints;
327 totpoints= numpoints+
dim+1;
328 else if (isregular) {
329 totpoints= numpoints;
332 totpoints += numpoints - 2;
333 }
else if (
dim == 3) {
335 totpoints += 2 * numpoints;
337 totpoints += 1 + numpoints;
342 totpoints= numpoints + isaxis;
343 totpoints += cubesize + diamondsize + addpoints;
344 totpoints += coincidentpoints*coincidenttotal;
348 for (s=command; *s; s++) {
349 if (issimplex2 && *s ==
'y')
355 }
else if (israndom) {
356 seed= (int)time(&timedata);
357 sprintf(seedbuf,
" t%d",
seed);
358 strncat(command, seedbuf,
sizeof(command)-strlen(command)-1);
359 t= strstr(command,
" t ");
369 NOcommand ?
"" : command,
371 qh->rbox_isinteger ?
"integer" :
"real");
379 if ((s= first_point)) {
387 if (isspace(*s) || !*s)
390 qh_fprintf_rbox(
qh,
qh->ferr, 6194,
"rbox error: missing comma after coordinate in %s\n\n", s);
395 for (k=
dim-count; k--; )
397 }
else if (count >
dim) {
398 qh_fprintf_rbox(
qh,
qh->ferr, 6195,
"rbox error: %d coordinates instead of %d coordinates in %s\n\n",
403 while ((s= strchr(s,
'P'))) {
411 if (issimplex+issimplex2) {
418 for (i=0; i<
dim; i++) {
419 for (k=0; k<
dim; k++)
420 *(simplexp++)= i==k ? 1.0 : 0.0;
422 for (k=0; k<
dim; k++)
425 for (i=0; i<
dim+1; i++) {
426 for (k=0; k<
dim; k++) {
428 *(simplexp++)= 2.0 * randr/randmax - 1.0;
434 for (i=0; i<
dim+1; i++) {
437 for (k=0; k<
dim; k++)
442 for (j=0; j<numpoints; j++) {
447 for (k=0; k<
dim; k++)
450 for (i=0; i<
dim+1; i++) {
452 factor= randr/randmax;
456 for (k=0; k<
dim; k++) {
457 simplexp= simplex + i*
dim + k;
458 coord[k] += factor * (*simplexp);
461 for (k=0; k<
dim; k++)
462 coord[k] *=
box/norm;
464 if(coincidentcount++ < coincidenttotal)
473 nthroot= (int)(pow((
double)numpoints, 1.0/
dim) + 0.99999);
476 for (i=0; i < numpoints; i++) {
478 for (k=0; k <
dim; k++) {
480 *(coordp++)= mult[0] * meshn + mult[1] * (-meshm);
482 *(coordp++)= mult[0] * meshm + mult[1] * meshn;
484 *(coordp++)= mult[k] * meshr;
487 if(coincidentcount++ < coincidenttotal)
489 for (k=0; k <
dim; k++) {
490 if (++mult[k] < nthroot)
497 else if (isregular && !islens) {
498 if (
dim != 2 &&
dim != 3) {
500 qh_fprintf_rbox(
qh,
qh->ferr, 6197,
"rbox error: regular points can be used only in 2-d and 3-d\n\n");
503 if (!isaxis || radius == 0.0) {
518 anglediff= 2.0 *
PI/numpoints;
519 for (i=0; i < numpoints; i++) {
521 x= radius * cos(angle);
522 y= radius * sin(angle);
528 norm= sqrt(1.0 +
x*
x +
y*
y);
535 norm= sqrt(1.0 +
x*
x +
y*
y);
544 else if (isregular && islens &&
dim == 2) {
548 anglediff= 2 * lensangle/(numpoints - 1);
549 cos_0= cos(lensangle);
550 for (i=0; i < numpoints; i++, angle -= anglediff) {
551 x= radius * sin(angle);
552 y= radius * (cos(angle) - cos_0);
556 if (i != 0 && i != numpoints - 1) {
564 else if (isregular && islens &&
dim != 2) {
567 qh_fprintf_rbox(
qh,
qh->ferr, 6198,
"rbox error: regular points can be used only in 2-d and 3-d\n\n");
571 anglediff= 2*
PI/numpoints;
576 offset= sqrt(radius * radius - (1-gap)*(1-gap)) - lensbase;
577 for (i=0; i < numpoints; i++, angle += anglediff) {
599 for (j=0; j < gendim; j++)
608 for (i=0; i < numpoints; i++) {
610 for (j=0; j < gendim; j++) {
612 coord[j]= 2.0 * randr/randmax - 1.0;
613 norm += coord[j] * coord[j];
623 rangap= 1.0 - gap * randr/randmax;
624 factor= radius * rangap / norm;
625 for (j=0; j<gendim; j++)
626 coord[j]= factor * coord[j];
628 }
else if (islens && issphere) {
634 rangap= 1.0 - gap * randr/randmax;
635 factor= rangap / norm;
636 for (j=0; j<gendim; j++)
637 coord[j]= factor * coord[j];
639 }
else if (islens && !issphere) {
646 coord[j]= -1.0 - coord[j] * gap;
648 coord[j]= 1.0 - coord[j] * gap;
650 }
else if (isspiral) {
653 qh_fprintf_rbox(
qh,
qh->ferr, 6199,
"rbox error: spiral distribution is available only in 3d\n\n");
656 coord[0]= cos(2*
PI*i/(numpoints - 1));
657 coord[1]= sin(2*
PI*i/(numpoints - 1));
658 coord[2]= 2.0*(double)i/(
double)(numpoints-1) - 1.0;
660 }
else if (issphere) {
664 factor *= 1.0 - width * randr/randmax;
666 for (j=0; j<
dim; j++)
667 coord[j]= factor * coord[j];
670 if (isaxis && issphere) {
673 for (j=0; j<gendim; j++)
674 norm += coord[j] * coord[j];
676 for (j=0; j<
dim; j++)
677 coord[j]= coord[j] / norm;
680 coord[
dim-1] *= 1 - width * randr/randmax;
683 }
else if (isaxis && !issphere) {
685 coord[
dim-1]= 2.0 * randr/randmax - 1.0;
688 coord[
dim-1]= lensbase;
689 for (j=0, norm= 0; j<
dim; j++)
690 norm += coord[j] * coord[j];
692 for (j=0; j<
dim; j++)
693 coord[j]= coord[j] * radius/ norm;
694 coord[
dim-1] -= lensbase;
697 coord[
dim-1] *= 1 - width * randr/randmax;
700 coord[
dim-1]= -coord[
dim-1];
702 }
else if (iswidth && !issphere) {
705 coord[j]= -1.0 - coord[j] * width;
707 coord[j]= 1.0 - coord[j] * width;
710 for (k=0; k<
dim; k++)
711 coord[k]= coord[k] *
box;
715 if(coincidentcount++ < coincidenttotal)
722 for (j=0; j<cubesize; j++) {
725 for (k=
dim-1; k>=0; k--) {
737 for (j=0; j<diamondsize; j++) {
740 for (k=
dim-1; k>=0; k--) {
765 if (
a - 0.5 < INT_MIN) {
766 qh_fprintf_rbox(
qh,
qh->ferr, 6200,
"rbox input error: negative coordinate %2.2g is too large. Reduce 'Bn'\n",
a);
769 return (
int)(
a - 0.5);
771 if (
a + 0.5 > INT_MAX) {
772 qh_fprintf_rbox(
qh,
qh->ferr, 6201,
"rbox input error: coordinate %2.2g is too large. Reduce 'Bn'\n",
a);
775 return (
int)(
a + 0.5);
781 if (
qh->rbox_isinteger)
789 if (
qh->rbox_isinteger)
797 if (
qh->rbox_isinteger)
809 for (k=0; k <
dim; k++)
820 for (i= 0; i<coincidentpoints; i++) {
824 for (k=0; k <
dim; k++) {
826 delta= 2.0 * randr/randmax - 1.0;
840 longjmp(
qh->rbox_errexit, exitcode);