← Back to branch summary

~ubuntu-branches/ubuntu/trusty/qhull/trusty-proposed

« back to all changes in this revision

Viewing changes to src/libqhull/rboxlib.c

  • Committer: Package Import Robot
  • Author(s): Barak A. Pearlmutter
  • Date: 2014-02-13 11:09:12 UTC
  • mfrom: (8.1.4 sid)
  • Revision ID: package-import@ubuntu.com-20140213110912-ifwyxorlsnnl1ebh
Tags: 2012.1-4
Add convenience link to #include <qhull/qhull.h> to simplify transition.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/libqhull/rboxlib.c
 
1
/*<html><pre>  -<a                             href="index.htm#TOC"
 
2
  >-------------------------------</a><a name="TOP">-</a>
 
3
 
 
4
   rboxlib.c
 
5
     Generate input points
 
6
 
 
7
   notes:
 
8
     For documentation, see prompt[] of rbox.c
 
9
     50 points generated for 'rbox D4'
 
10
 
 
11
   WARNING:
 
12
     incorrect range if qh_RANDOMmax is defined wrong (user.h)
 
13
*/
 
14
 
 
15
#include "random.h"
 
16
#include "libqhull.h"
 
17
 
 
18
#include <ctype.h>
 
19
#include <limits.h>
 
20
#include <math.h>
 
21
#include <setjmp.h>
 
22
#include <string.h>
 
23
#include <time.h>
 
24
#include <stdio.h>
 
25
#include <stdlib.h>
 
26
 
 
27
#ifdef _MSC_VER  /* Microsoft Visual C++ */
 
28
#pragma warning( disable : 4706)  /* assignment within conditional expression. */
 
29
#pragma warning( disable : 4996)  /* this function (strncat,sprintf,strcpy) or variable may be unsafe. */
 
30
#endif
 
31
 
 
32
#define MAXdim 200
 
33
#define PI 3.1415926535897932384
 
34
 
 
35
/* ------------------------------ prototypes ----------------*/
 
36
int roundi( double a);
 
37
void out1( double a);
 
38
void out2n( double a, double b);
 
39
void out3n( double a, double b, double c);
 
40
 
 
41
void    qh_fprintf_rbox(FILE *fp, int msgcode, const char *fmt, ... );
 
42
void    qh_free(void *mem);
 
43
void   *qh_malloc(size_t size);
 
44
int     qh_rand( void);
 
45
void    qh_srand( int seed);
 
46
 
 
47
 
 
48
/* ------------------------------ globals -------------------*/
 
49
 
 
50
/* No state is carried between rbox requests */
 
51
typedef struct rboxT rboxT;
 
52
struct rboxT {
 
53
  FILE *fout;
 
54
  FILE *ferr;
 
55
  int isinteger;
 
56
  double out_offset;
 
57
  jmp_buf errexit;        /* exit label for rboxpoints, defined by setjmp(), called by qh_errexit_rbox() */
 
58
};
 
59
 
 
60
 
 
61
int rbox_inuse= 0;
 
62
rboxT rbox;
 
63
 
 
64
/*-<a                             href="qh-qhull.htm#TOC"
 
65
  >-------------------------------</a><a name="rboxpoints">-</a>
 
66
 
 
67
  qh_rboxpoints( fout, ferr, rbox_command )
 
68
    Generate points to fout according to rbox options
 
69
    Report errors on ferr
 
70
 
 
71
  returns:
 
72
    0 (qh_ERRnone) on success
 
73
    1 (qh_ERRinput) on input error
 
74
    4 (qh_ERRmem) on memory error
 
75
    5 (qh_ERRqhull) on internal error
 
76
 
 
77
  notes:
 
78
    To avoid stdio, redefine qh_malloc, qh_free, and qh_fprintf_rbox (user.c)
 
79
    Rbox is not multithreaded.
 
80
 
 
81
  design:
 
82
    Straight line code (consider defining a struct and functions):
 
83
 
 
84
    Parse arguments into variables
 
85
    Determine the number of points
 
86
    Generate the points
 
87
*/
 
88
int qh_rboxpoints(FILE* fout, FILE* ferr, char* rbox_command) {
 
89
  int i,j,k;
 
90
  int gendim;
 
91
  int cubesize, diamondsize, seed=0, count, apex;
 
92
  int dim=3 , numpoints= 0, totpoints, addpoints=0;
 
93
  int issphere=0, isaxis=0,  iscdd= 0, islens= 0, isregular=0, iswidth=0, addcube=0;
 
94
  int isgap=0, isspiral=0, NOcommand= 0, adddiamond=0;
 
95
  int israndom=0, istime=0;
 
96
  int isbox=0, issimplex=0, issimplex2=0, ismesh=0;
 
97
  double width=0.0, gap=0.0, radius= 0.0;
 
98
  double coord[MAXdim], offset, meshm=3.0, meshn=4.0, meshr=5.0;
 
99
  double *simplex= NULL, *simplexp;
 
100
  int nthroot, mult[MAXdim];
 
101
  double norm, factor, randr, rangap, lensangle= 0, lensbase= 1;
 
102
  double anglediff, angle, x, y, cube= 0.0, diamond= 0.0;
 
103
  double box= qh_DEFAULTbox; /* scale all numbers before output */
 
104
  double randmax= qh_RANDOMmax;
 
105
  char command[200], seedbuf[200];
 
106
  char *s= command, *t, *first_point= NULL;
 
107
  time_t timedata;
 
108
  int exitcode;
 
109
 
 
110
  if (rbox_inuse) {
 
111
    qh_fprintf_rbox(rbox.ferr, 6188, "rbox error: rbox in use by another process.  Please lock calls to rbox.\n");
 
112
    return qh_ERRqhull;
 
113
  }
 
114
  rbox_inuse= True;
 
115
  rbox.ferr= ferr;
 
116
  rbox.fout= fout;
 
117
 
 
118
  exitcode= setjmp(rbox.errexit);
 
119
  if (exitcode) {
 
120
    /* same code for error exit and normal return */
 
121
    if (simplex)
 
122
        qh_free(simplex);
 
123
    rbox_inuse= False;
 
124
    return exitcode;
 
125
  }
 
126
 
 
127
  *command= '\0';
 
128
  strncat(command, rbox_command, sizeof(command)-strlen(command)-1);
 
129
 
 
130
  while (*s && !isspace(*s))  /* skip program name */
 
131
    s++;
 
132
  while (*s) {
 
133
    while (*s && isspace(*s))
 
134
      s++;
 
135
    if (*s == '-')
 
136
      s++;
 
137
    if (!*s)
 
138
      break;
 
139
    if (isdigit(*s)) {
 
140
      numpoints= qh_strtol(s, &s);
 
141
      continue;
 
142
    }
 
143
    /* ============= read flags =============== */
 
144
    switch (*s++) {
 
145
    case 'c':
 
146
      addcube= 1;
 
147
      t= s;
 
148
      while (isspace(*t))
 
149
        t++;
 
150
      if (*t == 'G')
 
151
        cube= qh_strtod(++t, &s);
 
152
      break;
 
153
    case 'd':
 
154
      adddiamond= 1;
 
155
      t= s;
 
156
      while (isspace(*t))
 
157
        t++;
 
158
      if (*t == 'G')
 
159
        diamond= qh_strtod(++t, &s);
 
160
      break;
 
161
    case 'h':
 
162
      iscdd= 1;
 
163
      break;
 
164
    case 'l':
 
165
      isspiral= 1;
 
166
      break;
 
167
    case 'n':
 
168
      NOcommand= 1;
 
169
      break;
 
170
    case 'r':
 
171
      isregular= 1;
 
172
      break;
 
173
    case 's':
 
174
      issphere= 1;
 
175
      break;
 
176
    case 't':
 
177
      istime= 1;
 
178
      if (isdigit(*s)) {
 
179
        seed= qh_strtol(s, &s);
 
180
        israndom= 0;
 
181
      }else
 
182
        israndom= 1;
 
183
      break;
 
184
    case 'x':
 
185
      issimplex= 1;
 
186
      break;
 
187
    case 'y':
 
188
      issimplex2= 1;
 
189
      break;
 
190
    case 'z':
 
191
      rbox.isinteger= 1;
 
192
      break;
 
193
    case 'B':
 
194
      box= qh_strtod(s, &s);
 
195
      isbox= 1;
 
196
      break;
 
197
    case 'D':
 
198
      dim= qh_strtol(s, &s);
 
199
      if (dim < 1
 
200
      || dim > MAXdim) {
 
201
        qh_fprintf_rbox(rbox.ferr, 6189, "rbox error: dimension, D%d, out of bounds (>=%d or <=0)", dim, MAXdim);
 
202
        qh_errexit_rbox(qh_ERRinput);
 
203
      }
 
204
      break;
 
205
    case 'G':
 
206
      if (isdigit(*s))
 
207
        gap= qh_strtod(s, &s);
 
208
      else
 
209
        gap= 0.5;
 
210
      isgap= 1;
 
211
      break;
 
212
    case 'L':
 
213
      if (isdigit(*s))
 
214
        radius= qh_strtod(s, &s);
 
215
      else
 
216
        radius= 10;
 
217
      islens= 1;
 
218
      break;
 
219
    case 'M':
 
220
      ismesh= 1;
 
221
      if (*s)
 
222
        meshn= qh_strtod(s, &s);
 
223
      if (*s == ',') {
 
224
        ++s;
 
225
        meshm= qh_strtod(s, &s);
 
226
      }else
 
227
        meshm= 0.0;
 
228
      if (*s == ',') {
 
229
        ++s;
 
230
        meshr= qh_strtod(s, &s);
 
231
      }else
 
232
        meshr= sqrt(meshn*meshn + meshm*meshm);
 
233
      if (*s && !isspace(*s)) {
 
234
        qh_fprintf_rbox(rbox.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;
 
236
      }
 
237
      break;
 
238
    case 'O':
 
239
      rbox.out_offset= qh_strtod(s, &s);
 
240
      break;
 
241
    case 'P':
 
242
      if (!first_point)
 
243
        first_point= s-1;
 
244
      addpoints++;
 
245
      while (*s && !isspace(*s))   /* read points later */
 
246
        s++;
 
247
      break;
 
248
    case 'W':
 
249
      width= qh_strtod(s, &s);
 
250
      iswidth= 1;
 
251
      break;
 
252
    case 'Z':
 
253
      if (isdigit(*s))
 
254
        radius= qh_strtod(s, &s);
 
255
      else
 
256
        radius= 1.0;
 
257
      isaxis= 1;
 
258
      break;
 
259
    default:
 
260
      qh_fprintf_rbox(rbox.ferr, 7070, "rbox error: unknown flag at %s.\nExecute 'rbox' without arguments for documentation.\n", s);
 
261
      qh_errexit_rbox(qh_ERRinput);
 
262
    }
 
263
    if (*s && !isspace(*s)) {
 
264
      qh_fprintf_rbox(rbox.ferr, 7071, "rbox error: missing space between flags at %s.\n", s);
 
265
      qh_errexit_rbox(qh_ERRinput);
 
266
    }
 
267
  }
 
268
 
 
269
  /* ============= defaults, constants, and sizes =============== */
 
270
  if (rbox.isinteger && !isbox)
 
271
    box= qh_DEFAULTzbox;
 
272
  if (addcube) {
 
273
    cubesize= (int)floor(ldexp(1.0,dim)+0.5);
 
274
    if (cube == 0.0)
 
275
      cube= box;
 
276
  }else
 
277
    cubesize= 0;
 
278
  if (adddiamond) {
 
279
    diamondsize= 2*dim;
 
280
    if (diamond == 0.0)
 
281
      diamond= box;
 
282
  }else
 
283
    diamondsize= 0;
 
284
  if (islens) {
 
285
    if (isaxis) {
 
286
        qh_fprintf_rbox(rbox.ferr, 6190, "rbox error: can not combine 'Ln' with 'Zn'\n");
 
287
        qh_errexit_rbox(qh_ERRinput);
 
288
    }
 
289
    if (radius <= 1.0) {
 
290
        qh_fprintf_rbox(rbox.ferr, 6191, "rbox error: lens radius %.2g should be greater than 1.0\n",
 
291
               radius);
 
292
        qh_errexit_rbox(qh_ERRinput);
 
293
    }
 
294
    lensangle= asin(1.0/radius);
 
295
    lensbase= radius * cos(lensangle);
 
296
  }
 
297
 
 
298
  if (!numpoints) {
 
299
    if (issimplex2)
 
300
        ; /* ok */
 
301
    else if (isregular + issimplex + islens + issphere + isaxis + isspiral + iswidth + ismesh) {
 
302
        qh_fprintf_rbox(rbox.ferr, 6192, "rbox error: missing count\n");
 
303
        qh_errexit_rbox(qh_ERRinput);
 
304
    }else if (adddiamond + addcube + addpoints)
 
305
        ; /* ok */
 
306
    else {
 
307
        numpoints= 50;  /* ./rbox D4 is the test case */
 
308
        issphere= 1;
 
309
    }
 
310
  }
 
311
  if ((issimplex + islens + isspiral + ismesh > 1)
 
312
  || (issimplex + issphere + isspiral + ismesh > 1)) {
 
313
    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");
 
314
    qh_errexit_rbox(qh_ERRinput);
 
315
  }
 
316
 
 
317
  /* ============= print header with total points =============== */
 
318
  if (issimplex || ismesh)
 
319
    totpoints= numpoints;
 
320
  else if (issimplex2)
 
321
    totpoints= numpoints+dim+1;
 
322
  else if (isregular) {
 
323
    totpoints= numpoints;
 
324
    if (dim == 2) {
 
325
        if (islens)
 
326
          totpoints += numpoints - 2;
 
327
    }else if (dim == 3) {
 
328
        if (islens)
 
329
          totpoints += 2 * numpoints;
 
330
      else if (isgap)
 
331
        totpoints += 1 + numpoints;
 
332
      else
 
333
        totpoints += 2;
 
334
    }
 
335
  }else
 
336
    totpoints= numpoints + isaxis;
 
337
  totpoints += cubesize + diamondsize + addpoints;
 
338
 
 
339
  /* ============= seed randoms =============== */
 
340
  if (istime == 0) {
 
341
    for (s=command; *s; s++) {
 
342
      if (issimplex2 && *s == 'y') /* make 'y' same seed as 'x' */
 
343
        i= 'x';
 
344
      else
 
345
        i= *s;
 
346
      seed= 11*seed + i;
 
347
    }
 
348
  }else if (israndom) {
 
349
    seed= (int)time(&timedata);
 
350
    sprintf(seedbuf, " t%d", seed);  /* appends an extra t, not worth removing */
 
351
    strncat(command, seedbuf, sizeof(command)-strlen(command)-1);
 
352
    t= strstr(command, " t ");
 
353
    if (t)
 
354
      strcpy(t+1, t+3); /* remove " t " */
 
355
  } /* else, seed explicitly set to n */
 
356
  qh_RANDOMseed_(seed);
 
357
 
 
358
  /* ============= print header =============== */
 
359
 
 
360
  if (iscdd)
 
361
      qh_fprintf_rbox(rbox.fout, 9391, "%s\nbegin\n        %d %d %s\n",
 
362
      NOcommand ? "" : command,
 
363
      totpoints, dim+1,
 
364
      rbox.isinteger ? "integer" : "real");
 
365
  else if (NOcommand)
 
366
      qh_fprintf_rbox(rbox.fout, 9392, "%d\n%d\n", dim, totpoints);
 
367
  else
 
368
      qh_fprintf_rbox(rbox.fout, 9393, "%d %s\n%d\n", dim, command, totpoints);
 
369
 
 
370
  /* ============= explicit points =============== */
 
371
  if ((s= first_point)) {
 
372
    while (s && *s) { /* 'P' */
 
373
      count= 0;
 
374
      if (iscdd)
 
375
        out1( 1.0);
 
376
      while (*++s) {
 
377
        out1( qh_strtod(s, &s));
 
378
        count++;
 
379
        if (isspace(*s) || !*s)
 
380
          break;
 
381
        if (*s != ',') {
 
382
          qh_fprintf_rbox(rbox.ferr, 6194, "rbox error: missing comma after coordinate in %s\n\n", s);
 
383
          qh_errexit_rbox(qh_ERRinput);
 
384
        }
 
385
      }
 
386
      if (count < dim) {
 
387
        for (k=dim-count; k--; )
 
388
          out1( 0.0);
 
389
      }else if (count > dim) {
 
390
        qh_fprintf_rbox(rbox.ferr, 6195, "rbox error: %d coordinates instead of %d coordinates in %s\n\n",
 
391
                  count, dim, s);
 
392
        qh_errexit_rbox(qh_ERRinput);
 
393
      }
 
394
      qh_fprintf_rbox(rbox.fout, 9394, "\n");
 
395
      while ((s= strchr(s, 'P'))) {
 
396
        if (isspace(s[-1]))
 
397
          break;
 
398
      }
 
399
    }
 
400
  }
 
401
 
 
402
  /* ============= simplex distribution =============== */
 
403
  if (issimplex+issimplex2) {
 
404
    if (!(simplex= (double*)qh_malloc( dim * (dim+1) * sizeof(double)))) {
 
405
      qh_fprintf_rbox(rbox.ferr, 6196, "rbox error: insufficient memory for simplex\n");
 
406
      qh_errexit_rbox(qh_ERRmem); /* qh_ERRmem */
 
407
    }
 
408
    simplexp= simplex;
 
409
    if (isregular) {
 
410
      for (i=0; i<dim; i++) {
 
411
        for (k=0; k<dim; k++)
 
412
          *(simplexp++)= i==k ? 1.0 : 0.0;
 
413
      }
 
414
      for (k=0; k<dim; k++)
 
415
        *(simplexp++)= -1.0;
 
416
    }else {
 
417
      for (i=0; i<dim+1; i++) {
 
418
        for (k=0; k<dim; k++) {
 
419
          randr= qh_RANDOMint;
 
420
          *(simplexp++)= 2.0 * randr/randmax - 1.0;
 
421
        }
 
422
      }
 
423
    }
 
424
    if (issimplex2) {
 
425
        simplexp= simplex;
 
426
      for (i=0; i<dim+1; i++) {
 
427
        if (iscdd)
 
428
          out1( 1.0);
 
429
        for (k=0; k<dim; k++)
 
430
          out1( *(simplexp++) * box);
 
431
        qh_fprintf_rbox(rbox.fout, 9395, "\n");
 
432
      }
 
433
    }
 
434
    for (j=0; j<numpoints; j++) {
 
435
      if (iswidth)
 
436
        apex= qh_RANDOMint % (dim+1);
 
437
      else
 
438
        apex= -1;
 
439
      for (k=0; k<dim; k++)
 
440
        coord[k]= 0.0;
 
441
      norm= 0.0;
 
442
      for (i=0; i<dim+1; i++) {
 
443
        randr= qh_RANDOMint;
 
444
        factor= randr/randmax;
 
445
        if (i == apex)
 
446
          factor *= width;
 
447
        norm += factor;
 
448
        for (k=0; k<dim; k++) {
 
449
          simplexp= simplex + i*dim + k;
 
450
          coord[k] += factor * (*simplexp);
 
451
        }
 
452
      }
 
453
      for (k=0; k<dim; k++)
 
454
        coord[k] /= norm;
 
455
      if (iscdd)
 
456
        out1( 1.0);
 
457
      for (k=0; k < dim; k++)
 
458
        out1( coord[k] * box);
 
459
      qh_fprintf_rbox(rbox.fout, 9396, "\n");
 
460
    }
 
461
    isregular= 0; /* continue with isbox */
 
462
    numpoints= 0;
 
463
  }
 
464
 
 
465
  /* ============= mesh distribution =============== */
 
466
  if (ismesh) {
 
467
    nthroot= (int)(pow((double)numpoints, 1.0/dim) + 0.99999);
 
468
    for (k=dim; k--; )
 
469
      mult[k]= 0;
 
470
    for (i=0; i < numpoints; i++) {
 
471
      for (k=0; k < dim; k++) {
 
472
        if (k == 0)
 
473
          out1( mult[0] * meshn + mult[1] * (-meshm));
 
474
        else if (k == 1)
 
475
          out1( mult[0] * meshm + mult[1] * meshn);
 
476
        else
 
477
          out1( mult[k] * meshr );
 
478
      }
 
479
      qh_fprintf_rbox(rbox.fout, 9397, "\n");
 
480
      for (k=0; k < dim; k++) {
 
481
        if (++mult[k] < nthroot)
 
482
          break;
 
483
        mult[k]= 0;
 
484
      }
 
485
    }
 
486
  }
 
487
  /* ============= regular points for 's' =============== */
 
488
  else if (isregular && !islens) {
 
489
    if (dim != 2 && dim != 3) {
 
490
      qh_fprintf_rbox(rbox.ferr, 6197, "rbox error: regular points can be used only in 2-d and 3-d\n\n");
 
491
      qh_errexit_rbox(qh_ERRinput);
 
492
    }
 
493
    if (!isaxis || radius == 0.0) {
 
494
      isaxis= 1;
 
495
      radius= 1.0;
 
496
    }
 
497
    if (dim == 3) {
 
498
      if (iscdd)
 
499
        out1( 1.0);
 
500
      out3n( 0.0, 0.0, -box);
 
501
      if (!isgap) {
 
502
        if (iscdd)
 
503
          out1( 1.0);
 
504
        out3n( 0.0, 0.0, box);
 
505
      }
 
506
    }
 
507
    angle= 0.0;
 
508
    anglediff= 2.0 * PI/numpoints;
 
509
    for (i=0; i < numpoints; i++) {
 
510
      angle += anglediff;
 
511
      x= radius * cos(angle);
 
512
      y= radius * sin(angle);
 
513
      if (dim == 2) {
 
514
        if (iscdd)
 
515
          out1( 1.0);
 
516
        out2n( x*box, y*box);
 
517
      }else {
 
518
        norm= sqrt(1.0 + x*x + y*y);
 
519
        if (iscdd)
 
520
          out1( 1.0);
 
521
        out3n( box*x/norm, box*y/norm, box/norm);
 
522
        if (isgap) {
 
523
          x *= 1-gap;
 
524
          y *= 1-gap;
 
525
          norm= sqrt(1.0 + x*x + y*y);
 
526
          if (iscdd)
 
527
            out1( 1.0);
 
528
          out3n( box*x/norm, box*y/norm, box/norm);
 
529
        }
 
530
      }
 
531
    }
 
532
  }
 
533
  /* ============= regular points for 'r Ln D2' =============== */
 
534
  else if (isregular && islens && dim == 2) {
 
535
    double cos_0;
 
536
 
 
537
    angle= lensangle;
 
538
    anglediff= 2 * lensangle/(numpoints - 1);
 
539
    cos_0= cos(lensangle);
 
540
    for (i=0; i < numpoints; i++, angle -= anglediff) {
 
541
      x= radius * sin(angle);
 
542
      y= radius * (cos(angle) - cos_0);
 
543
      if (iscdd)
 
544
        out1( 1.0);
 
545
      out2n( x*box, y*box);
 
546
      if (i != 0 && i != numpoints - 1) {
 
547
        if (iscdd)
 
548
          out1( 1.0);
 
549
        out2n( x*box, -y*box);
 
550
      }
 
551
    }
 
552
  }
 
553
  /* ============= regular points for 'r Ln D3' =============== */
 
554
  else if (isregular && islens && dim != 2) {
 
555
    if (dim != 3) {
 
556
      qh_fprintf_rbox(rbox.ferr, 6198, "rbox error: regular points can be used only in 2-d and 3-d\n\n");
 
557
      qh_errexit_rbox(qh_ERRinput);
 
558
    }
 
559
    angle= 0.0;
 
560
    anglediff= 2* PI/numpoints;
 
561
    if (!isgap) {
 
562
      isgap= 1;
 
563
      gap= 0.5;
 
564
    }
 
565
    offset= sqrt(radius * radius - (1-gap)*(1-gap)) - lensbase;
 
566
    for (i=0; i < numpoints; i++, angle += anglediff) {
 
567
      x= cos(angle);
 
568
      y= sin(angle);
 
569
      if (iscdd)
 
570
        out1( 1.0);
 
571
      out3n( box*x, box*y, 0.0);
 
572
      x *= 1-gap;
 
573
      y *= 1-gap;
 
574
      if (iscdd)
 
575
        out1( 1.0);
 
576
      out3n( box*x, box*y, box * offset);
 
577
      if (iscdd)
 
578
        out1( 1.0);
 
579
      out3n( box*x, box*y, -box * offset);
 
580
    }
 
581
  }
 
582
  /* ============= apex of 'Zn' distribution + gendim =============== */
 
583
  else {
 
584
    if (isaxis) {
 
585
      gendim= dim-1;
 
586
      if (iscdd)
 
587
        out1( 1.0);
 
588
      for (j=0; j < gendim; j++)
 
589
        out1( 0.0);
 
590
      out1( -box);
 
591
      qh_fprintf_rbox(rbox.fout, 9398, "\n");
 
592
    }else if (islens)
 
593
      gendim= dim-1;
 
594
    else
 
595
      gendim= dim;
 
596
    /* ============= generate random point in unit cube =============== */
 
597
    for (i=0; i < numpoints; i++) {
 
598
      norm= 0.0;
 
599
      for (j=0; j < gendim; j++) {
 
600
        randr= qh_RANDOMint;
 
601
        coord[j]= 2.0 * randr/randmax - 1.0;
 
602
        norm += coord[j] * coord[j];
 
603
      }
 
604
      norm= sqrt(norm);
 
605
      /* ============= dim-1 point of 'Zn' distribution ========== */
 
606
      if (isaxis) {
 
607
        if (!isgap) {
 
608
          isgap= 1;
 
609
          gap= 1.0;
 
610
        }
 
611
        randr= qh_RANDOMint;
 
612
        rangap= 1.0 - gap * randr/randmax;
 
613
        factor= radius * rangap / norm;
 
614
        for (j=0; j<gendim; j++)
 
615
          coord[j]= factor * coord[j];
 
616
      /* ============= dim-1 point of 'Ln s' distribution =========== */
 
617
      }else if (islens && issphere) {
 
618
        if (!isgap) {
 
619
          isgap= 1;
 
620
          gap= 1.0;
 
621
        }
 
622
        randr= qh_RANDOMint;
 
623
        rangap= 1.0 - gap * randr/randmax;
 
624
        factor= rangap / norm;
 
625
        for (j=0; j<gendim; j++)
 
626
          coord[j]= factor * coord[j];
 
627
      /* ============= dim-1 point of 'Ln' distribution ========== */
 
628
      }else if (islens && !issphere) {
 
629
        if (!isgap) {
 
630
          isgap= 1;
 
631
          gap= 1.0;
 
632
        }
 
633
        j= qh_RANDOMint % gendim;
 
634
        if (coord[j] < 0)
 
635
          coord[j]= -1.0 - coord[j] * gap;
 
636
        else
 
637
          coord[j]= 1.0 - coord[j] * gap;
 
638
      /* ============= point of 'l' distribution =============== */
 
639
      }else if (isspiral) {
 
640
        if (dim != 3) {
 
641
          qh_fprintf_rbox(rbox.ferr, 6199, "rbox error: spiral distribution is available only in 3d\n\n");
 
642
          longjmp(rbox.errexit,qh_ERRinput);
 
643
        }
 
644
        coord[0]= cos(2*PI*i/(numpoints - 1));
 
645
        coord[1]= sin(2*PI*i/(numpoints - 1));
 
646
        coord[2]= 2.0*(double)i/(double)(numpoints-1) - 1.0;
 
647
      /* ============= point of 's' distribution =============== */
 
648
      }else if (issphere) {
 
649
        factor= 1.0/norm;
 
650
        if (iswidth) {
 
651
          randr= qh_RANDOMint;
 
652
          factor *= 1.0 - width * randr/randmax;
 
653
        }
 
654
        for (j=0; j<dim; j++)
 
655
          coord[j]= factor * coord[j];
 
656
      }
 
657
      /* ============= project 'Zn s' point in to sphere =============== */
 
658
      if (isaxis && issphere) {
 
659
        coord[dim-1]= 1.0;
 
660
        norm= 1.0;
 
661
        for (j=0; j<gendim; j++)
 
662
          norm += coord[j] * coord[j];
 
663
        norm= sqrt(norm);
 
664
        for (j=0; j<dim; j++)
 
665
          coord[j]= coord[j] / norm;
 
666
        if (iswidth) {
 
667
          randr= qh_RANDOMint;
 
668
          coord[dim-1] *= 1 - width * randr/randmax;
 
669
        }
 
670
      /* ============= project 'Zn' point onto cube =============== */
 
671
      }else if (isaxis && !issphere) {  /* not very interesting */
 
672
        randr= qh_RANDOMint;
 
673
        coord[dim-1]= 2.0 * randr/randmax - 1.0;
 
674
      /* ============= project 'Ln' point out to sphere =============== */
 
675
      }else if (islens) {
 
676
        coord[dim-1]= lensbase;
 
677
        for (j=0, norm= 0; j<dim; j++)
 
678
          norm += coord[j] * coord[j];
 
679
        norm= sqrt(norm);
 
680
        for (j=0; j<dim; j++)
 
681
          coord[j]= coord[j] * radius/ norm;
 
682
        coord[dim-1] -= lensbase;
 
683
        if (iswidth) {
 
684
          randr= qh_RANDOMint;
 
685
          coord[dim-1] *= 1 - width * randr/randmax;
 
686
        }
 
687
        if (qh_RANDOMint > randmax/2)
 
688
          coord[dim-1]= -coord[dim-1];
 
689
      /* ============= project 'Wn' point toward boundary =============== */
 
690
      }else if (iswidth && !issphere) {
 
691
        j= qh_RANDOMint % gendim;
 
692
        if (coord[j] < 0)
 
693
          coord[j]= -1.0 - coord[j] * width;
 
694
        else
 
695
          coord[j]= 1.0 - coord[j] * width;
 
696
      }
 
697
      /* ============= write point =============== */
 
698
      if (iscdd)
 
699
        out1( 1.0);
 
700
      for (k=0; k < dim; k++)
 
701
        out1( coord[k] * box);
 
702
      qh_fprintf_rbox(rbox.fout, 9399, "\n");
 
703
    }
 
704
  }
 
705
 
 
706
  /* ============= write cube vertices =============== */
 
707
  if (addcube) {
 
708
    for (j=0; j<cubesize; j++) {
 
709
      if (iscdd)
 
710
        out1( 1.0);
 
711
      for (k=dim-1; k>=0; k--) {
 
712
        if (j & ( 1 << k))
 
713
          out1( cube);
 
714
        else
 
715
          out1( -cube);
 
716
      }
 
717
      qh_fprintf_rbox(rbox.fout, 9400, "\n");
 
718
    }
 
719
  }
 
720
 
 
721
  /* ============= write diamond vertices =============== */
 
722
  if (adddiamond) {
 
723
    for (j=0; j<diamondsize; j++) {
 
724
      if (iscdd)
 
725
        out1( 1.0);
 
726
      for (k=dim-1; k>=0; k--) {
 
727
        if (j/2 != k)
 
728
          out1( 0.0);
 
729
        else if (j & 0x1)
 
730
          out1( diamond);
 
731
        else
 
732
          out1( -diamond);
 
733
      }
 
734
      qh_fprintf_rbox(rbox.fout, 9401, "\n");
 
735
    }
 
736
  }
 
737
 
 
738
  if (iscdd)
 
739
    qh_fprintf_rbox(rbox.fout, 9402, "end\nhull\n");
 
740
 
 
741
  /* same code for error exit and normal return */
 
742
  if (simplex)
 
743
    qh_free(simplex);
 
744
  rbox_inuse= False;
 
745
  return qh_ERRnone;
 
746
} /* rboxpoints */
 
747
 
 
748
/*------------------------------------------------
 
749
outxxx - output functions
 
750
*/
 
751
int roundi( double a) {
 
752
  if (a < 0.0) {
 
753
    if (a - 0.5 < INT_MIN) {
 
754
      qh_fprintf_rbox(rbox.ferr, 6200, "rbox input error: negative coordinate %2.2g is too large.  Reduce 'Bn'\n", a);
 
755
      qh_errexit_rbox(qh_ERRinput);
 
756
    }
 
757
    return (int)(a - 0.5);
 
758
  }else {
 
759
    if (a + 0.5 > INT_MAX) {
 
760
      qh_fprintf_rbox(rbox.ferr, 6201, "rbox input error: coordinate %2.2g is too large.  Reduce 'Bn'\n", a);
 
761
      qh_errexit_rbox(qh_ERRinput);
 
762
    }
 
763
    return (int)(a + 0.5);
 
764
  }
 
765
} /* roundi */
 
766
 
 
767
void out1(double a) {
 
768
 
 
769
  if (rbox.isinteger)
 
770
    qh_fprintf_rbox(rbox.fout, 9403, "%d ", roundi( a+rbox.out_offset));
 
771
  else
 
772
    qh_fprintf_rbox(rbox.fout, 9404, qh_REAL_1, a+rbox.out_offset);
 
773
} /* out1 */
 
774
 
 
775
void out2n( double a, double b) {
 
776
 
 
777
  if (rbox.isinteger)
 
778
    qh_fprintf_rbox(rbox.fout, 9405, "%d %d\n", roundi(a+rbox.out_offset), roundi(b+rbox.out_offset));
 
779
  else
 
780
    qh_fprintf_rbox(rbox.fout, 9406, qh_REAL_2n, a+rbox.out_offset, b+rbox.out_offset);
 
781
} /* out2n */
 
782
 
 
783
void out3n( double a, double b, double c) {
 
784
 
 
785
  if (rbox.isinteger)
 
786
    qh_fprintf_rbox(rbox.fout, 9407, "%d %d %d\n", roundi(a+rbox.out_offset), roundi(b+rbox.out_offset), roundi(c+rbox.out_offset));
 
787
  else
 
788
    qh_fprintf_rbox(rbox.fout, 9408, qh_REAL_3n, a+rbox.out_offset, b+rbox.out_offset, c+rbox.out_offset);
 
789
} /* out3n */
 
790
 
 
791
void qh_errexit_rbox(int exitcode)
 
792
{
 
793
    longjmp(rbox.errexit, exitcode);
 
794
} /* rbox_errexit */
 
795