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);