/* ***** PLANES.C ***** */ #include #include #include #include /* This is a program I wrote to generate the paper */ /* "cut-out" that folds up into a 3-D polytope. */ /* The input is "half-space cutting planes" and */ /* the output is an input file for the DM program */ /* which sees "L K X1 Y1 X2 Y2" as a line of color */ /* K from (X1,Y1) to (X2,Y2). Lines beginning with */ /* a pound sign (#) are comments to Draftsman-EE. */ /* I've added a section for personal computers with */ /* an EGA display when compiled with PC set to one. */ /* # This is the input file for my example. */ /* # */ /* A -1 0 0 0 */ /* B 0 -1 0 0 */ /* C 0 0 -1 0 */ /* D 55 -13 77 539 */ /* E -7 2 -14 14 */ /* F 8 2 4 64 */ /* G 3 9 -8 78 */ /* H -24 12 36 252 */ /* I 1 -2 -1 9 */ /* J -3 6 -2 42 */ /* K -4 -2 -1 0 */ /* L 9 3 1 75 */ /* This program follows the "Adam" style closely. */ /* Maximum line-length is 72 with a tab stop of 8. */ /* Square brackets count as three characters. */ /* Anything after a tab on a line is a comment. */ /* Code should be written for EBCDIC as well as */ /* ASCII where practical; use ISUPPER, for example, */ /* instead of comparing to 'A' and 'Z' constants. */ /* The words "err" "or" and "war" "ning" and */ /* slash-star and star-slash not as comments do */ /* not appear as one contiguous string. You can */ /* break them up in character strings as above */ /* and use "err*r" and "warn*ng" in comments. */ /* The restriction applies to upper and lower case. */ /* Indent is two characters for #IF and braces. */ /* Indenting should communicate program intent */ /* even if it breaks other "reasonable" rules. */ /* Braces go on their own line (unless they fit on */ /* one line) and indent less than their contents. */ /* It's considered good to indicate where matching */ /* braces can be found; I usually use line numbers. */ /* Labels (as in GOTO statements) also are marked. */ /* All functions and subroutines are prototyped. */ /* Source files, functions, and subroutines all */ /* have a specific 2-tab five-star header comment */ /* where 16 spaces replace the two tabs in EBCDIC. */ /* (I have tools that use those to find functions.) */ /* Also, all variables are declared at the head of */ /* a source file (or included header files there) */ /* or at the start of functions and subroutines. */ /* I don't like surprises declared mid-stream. */ /* Memory allocation is done once (and not undone) */ /* and is in its own separate subroutine. The only */ /* data structure is the array and all arrays that */ /* are not character strings have bounds checking. */ /* The subscript check can be against the number of */ /* items actually in use rather than maximum size. */ /* I'm not allergic to the comma operator or the */ /* GOTO statement. Labels and TYPEDEFs are first- */ /* letter-of-each-word upper case, variables and */ /* array macros are lower case, other #DEFINEs and */ /* ENUMs are all upper case. Two things can be the */ /* same except for case if they refer to the exact */ /* same thing; array macros are my obvious example. */ /* Variables in comments are written in upper case. */ /* Comma operands should be on the same line or */ /* visually obvious via indenting. */ /* Subscripts of arrays or array macros get no */ /* space on either side of commas or parentheses. */ /* Assignments and comparisons get space on both */ /* sides and parentheses get space on the outside */ /* (and IF and WHILE get at least two on the right) */ /* except for the left paren of an array or macro. */ /* Commas get space on the right most of the time. */ /* Bit operators (&, |, ^) and Booleans (&&, ||) */ /* get space on both sides and parentheses are */ /* used to define precedence even when unnecessary. */ /* I'm in favor of extra parentheses for clarity. */ /* I don't do assignments inside the decision part */ /* of IF, FOR, or WHILE statements and I stay away */ /* from testing err*r conditions by putting SCANF */ /* or any other subroutine in the decision part. */ /* A TRUE/FALSE function, however, belongs there. */ /* (If you *must* assign something in a decision */ /* area, then it should be inside extra parens.) */ /* I emphatically avoid the use of ! for NOT as it */ /* is often overlooked and I don't like "IF (J)" */ /* for "IF (J != 0)" unless J is TRUE/FALSE. */ /* Increment (++) and decrement (--) operators are */ /* separate from other operations so there is no */ /* pre- or post- increment or decrement issue. */ /* Also, I prefer P[0] to *P for the first char. */ /* First letters describe variable type as we used */ /* to do in FORTRAN. P is a pointer, O, Q and Z */ /* are character strings but Q can be a character, */ /* I through N are integers, and almost anything */ /* else is some kind of floating point. */ /* FILE pointers begin with X and sometimes Y. */ /* TRUE/FALSE, shorts, and longs are all integers. */ /* Single and double precision are floating point. */ /* Array pointers are named for what they point to. */ /* And you should seldom need to declare much else! */ /* Each primitive gets its own letter (or number). */ /* Most primitives have arrays but some don't. */ /* Scalar primitive variables have primitive */ /* letter second, so JP is (integer) plane, KV is */ /* (integer) vertex, FX would be (float) x-coef. */ /* Array names end in O (for "of") and a list of */ /* primitives which the subscripts should match. */ /* The number of primitive P in use is JNOP. */ /* The "memory max" size of primitive P is MEMP. */ /* The dimensioned size of primitive P is NUMP. */ /* (I have tools that use this to allocate memory.) */ /* Single letter variables can be anything. */ /* When obvious, so can double letter variables. */ /* Variable names longer than eight are TRUE/FALSE. */ /* Only subroutines have names longer than eight. */ /* (A C function becomes a "subroutine" to me when */ /* it "returns a void" or the return is some kind */ /* of status conditions, so PRINTF is a subroutine. */ /* List of primitives. */ #define NUML 210 /* The letter L stands for output lines. */ #define NUMP 30 /* The letter P stands for planes. */ #define NUMV 100 /* The letter V stands for vertices. */ /* The letter X stands for x-coordinate. */ /* The letter Y stands for y-coordinate. */ /* The letter Z stands for z-coordinate. */ #define FALSE 0 #define TRUE 1 #define SAME 0 #define EPS 0.000001 #define NOWHERE -123456.78 /* Array bounds checking. In C there is no */ /* "automatic" range checker like the ones that */ /* usually come with FORTRAN or PASCAL compilers. */ /* All arrays except for character strings should */ /* have bounds checking able to be turned on. */ #ifndef RANGE #define RANGE 0 #endif #if RANGE #define Range(n,nmin,nmax) zrange(n,nmin,nmax,__FILE__,__LINE__,#n) #else #define Range(x,y,z) (x) #endif /* This variable turns on EGA graphics stuff. */ #ifndef PC #define PC 0 #endif #if PC #include #include /* For 80386 computers, use INT386 instead of INT86. */ #ifdef __386__ #define int86 int386 #endif #endif /* * * * * * Lines * * * * * */ #define fxol(jl) Fxol[Range(jl,1,jnol)] double *Fxol; /* x-coordinate of first endpoint */ #define fyol(jl) Fyol[Range(jl,1,jnol)] double *Fyol; /* y-coordinate of first endpoint */ #define gxol(jl) Gxol[Range(jl,1,jnol)] double *Gxol; /* x-coordinate of second endpoint */ #define gyol(jl) Gyol[Range(jl,1,jnol)] double *Gyol; /* y-coordinate of second endpoint */ #define npol(jl) Npol[Range(jl,1,jnol)] int *Npol; /* plane line is associated with */ int jnol; /* number of lines */ /* * * * * * Planes * * * * * */ #define qnop(jp) Qnop[Range(jp,1,jnop)] char *Qnop; /* name of plane */ #define fxop(jp) Fxop[Range(jp,1,jnop)] double *Fxop; /* x coefficient of plane */ #define fyop(jp) Fyop[Range(jp,1,jnop)] double *Fyop; /* y coefficient of plane */ #define fzop(jp) Fzop[Range(jp,1,jnop)] double *Fzop; /* z coefficient of plane */ #define fbop(jp) Fbop[Range(jp,1,jnop)] double *Fbop; /* right hand side */ #define iop(jp) Iop[Range(jp,1,jnop)] int *Iop; /* flag variable */ int jnop; /* number of planes */ /* * * * * * Vertices * * * * * */ #define fxov(jv) Fxov[Range(jv,1,jnov)] double *Fxov; /* x coefficient of vertex */ #define fyov(jv) Fyov[Range(jv,1,jnov)] double *Fyov; /* y coefficient of vertex */ #define fzov(jv) Fzov[Range(jv,1,jnov)] double *Fzov; /* z coefficient of vertex */ #define gxov(jv) Gxov[Range(jv,1,jnov)] double *Gxov; /* x coefficient in 2-D */ #define gyov(jv) Gyov[Range(jv,1,jnov)] double *Gyov; /* y coefficient in 2-D */ #define hxov(jv) Hxov[Range(jv,1,jnov)] double *Hxov; /* x coef. in foldout */ #define hyov(jv) Hyov[Range(jv,1,jnov)] double *Hyov; /* y coef. in foldout */ int jnov; /* number of vertices */ /* * * * * * Plane by Vertex * * * * * */ #define nopv(jp,jv) Nopv[Range(jp,1,jnop)][Range(jv,1,jnov)] int **Nopv; /* one if V in P */ #define mopv(jp,jv) Mopv[Range(jp,1,jnop)][Range(jv,1,jnov)] int **Mopv; /* 1 if V in P, 2D */ /* * * * * * Global variables * * * * * */ int idebug; double det; double a, b, c, d, e, f, g, h, o; double aa, bb, cc, dd, ee, ff, gg, hh, oo; double r, s, t, u, v, w, x, y, z; /* * * * * * Subroutines and functions * * * * * */ void draw_picture (void); int iclose (double a, double b); void invert (void); void *mal (int n); void memory_allocation (void); void watch (void); int zrange (int m, int min, int max, char *ofile, int line, char *p); void zzrange (int i); /* ***** MAIN ***** */ int main (int argc, char **argv) { /* 293-709 */ char *p, q[256]; /* generic character string */ int j; /* generic integer variable */ int jp, kp, lp, mp; /* planes */ int jv, kv, lv, mv; /* vertices */ int i_did_something; FILE *xfile, *yfile; /* input and output file pointers */ /* Check command line for /DEBUG option at end. */ p = argv[argc-1]; if (strcmp (p, "/debug") == SAME) idebug = TRUE, argc--; else idebug = FALSE; /* Find input file from command line or terminal. */ if (argc < 2) printf (" File name for input --> "), gets (q); else strcpy (q, argv[1]); xfile = fopen (q, "r"); if (xfile != 0) printf ("Opened file %s.\n", q); else printf ("Cannot open file %s.\n", q), exit (9); /* Find output file from command line or terminal. */ if (argc < 3) printf ("File name for output --> "), gets (q); else strcpy (q, argv[2]); yfile = fopen (q, "w"); if (yfile != 0) printf ("Opened file %s.\n", q); else printf ("Cannot open file %s.\n", q), exit (9); /* Allocate array memory. */ memory_allocation (); /* The reason most of my output lines start with */ /* a pound sign (#) is that the DM-EE program I */ /* use to print the lines considers any line with */ /* a pound sign as a comment and only the lines */ /* starting with L (below) are part of the display. */ /* Read in input planes. */ jnop = 0; while (TRUE) { /* 337-362 */ fgets (q, 256, xfile); /* read line from file */ if (feof (xfile)) break; /* if end-of-file, then stop reading */ if (q[0] == '#') continue; if ((q[0] == 'e') && (q[1] == 'n')) break; jnop++; if (jnop > NUMP) { /* 346-349 */ printf ("Too many planes (%d)\n", jnop); exit (9); } /* 346-349 */ j = sscanf (q, "%c %lf %lf %lf %lf", &qnop(jnop), &fxop(jnop), &fyop(jnop), &fzop(jnop), &fbop(jnop)); if (j != 5) { /* 354-357 */ printf ("Format err" "or in %s.\n", q); exit (9); } /* 354-357 */ sprintf (q, "# Plane=%c %8.2lfx +%8.2lfy +%8.2lfz <= %8.2lf", qnop(jnop), fxop(jnop), fyop(jnop), fzop(jnop), fbop(jnop)); printf ("%s\n", q), fprintf (yfile, "%s\n", q); } /* 337-362 */ fclose (xfile); /* Tell user what we read. */ printf ("%d planes read from input file.\n", jnop); fprintf (yfile, "#\n"); /* Now find all vertices. */ jnov = 0; for (jp = 1; jp <= jnop-2; jp++) { /* 373-433 */ for (kp = jp+1; kp <= jnop-1; kp++) { /* 375-432 */ for (lp = kp+1; lp <= jnop; lp++) { /* 377-431 */ /* Find matrix A=(a,b,c,d,e,f,g,h,o). */ a = fxop(jp), b = fyop(jp), c = fzop(jp); d = fxop(kp), e = fyop(kp), f = fzop(kp); g = fxop(lp), h = fyop(lp), o = fzop(lp); /* Find right hand side vector (r,s,t). */ r = fbop(jp), s = fbop(kp), t = fbop(lp); /* Try to invert matrix A. */ invert (); if (det == 0.0) continue; /* no inverse, no vertex */ /* Multiple A-inverse by right hand side. */ x = (aa*r)+(bb*s)+(cc*t); y = (dd*r)+(ee*s)+(ff*t); z = (gg*r)+(hh*s)+(oo*t); /* Is it feasible? */ for (mp = 1; mp <= jnop; mp++) { /* 397-400 */ u = (fxop(mp)*x)+(fyop(mp)*y)+(fzop(mp)*z); if (u > fbop(mp)+EPS) break; } /* 397-400 */ if (mp <= jnop) continue; /* Not feasible. */ /* Is this already a vertex? */ for (jv = 1; jv <= jnov; jv++) { /* 405-409 */ if ((iclose (fxov(jv), x)) && (iclose (fyov(jv), y)) && (iclose (fzov(jv), z))) break; } /* 405-409 */ if (jv > jnov) /* No, add a new vertex. */ { /* 412-424 */ jnov++; if (jnov > NUMV) { /* 415-418 */ printf ("Too many vertices (%d)", jnov); exit (9); } /* 415-418 */ jv = jnov; fxov(jv) = x, fyov(jv) = y, fzov(jv) = z; for (mp = 1; mp <= jnop; mp++) nopv(mp,jv) = mopv(mp,jv) = FALSE; } /* 412-424 */ /* The vertex JV is in all three planes. */ nopv(jp,jv) = TRUE; nopv(kp,jv) = TRUE; nopv(lp,jv) = TRUE; } /* 377-431 */ } /* 375-432 */ } /* 373-433 */ /* Print the vertices and their planes. */ for (jv = 1; jv <= jnov; jv++) { /* 438-447 */ j = sprintf (q, "# Vertex %d = %8.2lf %8.2lf %8.2lf ", jv, fxov(jv), fyov(jv), fzov(jv)); for (mp = 1; mp <= jnop; mp++) { /* 442-445 */ if (nopv(mp,jv)) j += sprintf (q+j, " %c", qnop(mp)); else j += sprintf (q+j, " "); } /* 442-445 */ printf ("%s\n", q), fprintf (yfile, "%s\n", q); } /* 438-447 */ fprintf (yfile, "#\n"); /* Now for each plane, we figure out a two */ /* dimensional polygon "cutout" and then fit */ /* it into a foldout to print in DMEE format. */ /* We haven't done any planes yet. */ for (jp = 1; jp <= jnop; jp++) iop(jp) = FALSE; /* All vertices start nowhere at all. */ for (jv = 1; jv <= jnov; jv++) gxov(jv) = gyov(jv) = hxov(jv) = hyov(jv) = NOWHERE; jnol = 0, i_did_something = TRUE; while (i_did_something) { /* 464-702 */ i_did_something = FALSE; for (jp = 1; jp <= jnop; jp++) { /* 468-701 */ /* Have we already done this plane? */ if (iop(jp)) continue; /* The first plane is always OK to do. */ /* Otherwise, at least two vertices of this plane */ /* must be in a plane processed already. */ if (jp != 1) { /* 475-487 */ for (kp = 1; kp <= jnop; kp++) { /* 477-485 */ if (iop(kp) == FALSE) continue; j = 0; for (jv = 1; jv <= jnov; jv++) if ((nopv(jp,jv)) && (mopv(kp,jv))) j++, lv = kv, kv = jv; if (j >= 2) break; } /* 477-485 */ if (kp > jnop) continue; } /* 475-487 */ iop(jp) = i_did_something = TRUE; /* We need an orthogonal matrix representing */ /* a basis where the first coefficient is constant */ /* for all points in plane JP, including vertices. */ /* We form three vectors, U, V, and W. */ /* U=(a,b,c) is the plane's normal direction. */ /* V and W are handy perpendicular vectors. */ a = fxop(jp), b = fyop(jp), c = fzop(jp); r = sqrt ((a*a)+(b*b)+(c*c)); a /= r, b /= r, c /= r; /* That's one. */ if ((a > -0.7) && (a < 0.7)) d = 1.0, e = 0.0, f = 0.0; else d = 0.0, e = 1.0, f = 0.0; s = (a*d)+(b*e)+(c*f); /* Graham-Schmidt projection */ d -= s*a, e -= s*b, f -= s*c; /* Graham-Schmidt subtraction */ r = sqrt ((d*d)+(e*e)+(f*f)); /* Find length to normalize. */ d /= r, e /= r, f /= r; /* That's two. */ /* Use the cross product to preserve orientation. */ g = (b*f)-(c*e), h = (c*d)-(a*f), o = (a*e)-(b*d); s = (a*g)+(b*h)+(c*o); /* G-S projection, one */ g -= s*a, h -= s*b, o -= s*c; /* G-S subtraction, one */ s = (d*g)+(e*h)+(f*o); /* G-S projection, two */ g -= s*d, h -= s*e, o -= s*f; /* G-S subtraction, two */ r = sqrt ((g*g)+(h*h)+(o*o)); /* Find length to normalize. */ g /= r, h /= r, o /= r; /* That's all three. */ /* Print results if DEBUG flag is turned on. */ if (idebug) { /* 523-536 */ sprintf (q, "#\n"); printf ("%s", q), fprintf (yfile, "%s", q); sprintf (q, "# Plane %c.\n", qnop(jp)); printf ("%s", q), fprintf (yfile, "%s", q); sprintf (q, "#\n"); printf ("%s", q), fprintf (yfile, "%s", q); sprintf (q, "# a b c %8.2lf %8.2lf %8.2lf\n", a, b, c); printf ("%s", q), fprintf (yfile, "%s", q); sprintf (q, "# d e f %8.2lf %8.2lf %8.2lf\n", d, e, f); printf ("%s", q), fprintf (yfile, "%s", q); sprintf (q, "# g h o %8.2lf %8.2lf %8.2lf\n", g, h, o); printf ("%s", q), fprintf (yfile, "%s", q); } /* 523-536 */ for (jv = 1; jv <= jnov; jv++) { /* 539-558 */ /* We only want vertices in plane JP. */ if (nopv(jp,jv) == FALSE) continue; /* Find vertex vector (R,S,T). */ r = fxov(jv), s = fyov(jv), t = fzov(jv); /* Transform by above matrix. */ u = (a*r)+(b*s)+(c*t); v = gxov(jv) = (d*r)+(e*s)+(f*t); w = gyov(jv) = (g*r)+(h*s)+(o*t); /* Print results if DEBUG flag is turned on. */ if (idebug) { /* 550-557 */ sprintf (q, "#%d\n", jv); printf ("%s", q), fprintf (yfile, "%s", q); sprintf (q, "# r s t %8.2lf %8.2lf %8.2lf\n", r, s, t); printf ("%s", q), fprintf (yfile, "%s", q); sprintf (q, "# u v w %8.2lf %8.2lf %8.2lf\n", u, v, w); printf ("%s", q), fprintf (yfile, "%s", q); } /* 550-557 */ } /* 539-558 */ /* The first plane maps 2-D to folding identically. */ if (jp == 1) { /* 562-581 */ for (jv = 1; jv <= jnov; jv++) { /* 564-579 */ if (nopv(jp,jv)) { /* 566-578 */ x = hxov(jv) = gxov(jv); y = hyov(jv) = gyov(jv); nopv(jp,jv) = mopv(jp,jv) = TRUE; if (idebug) { /* 572-577 */ sprintf (q, "#\n"); printf ("%s", q), fprintf (yfile, "%s", q); sprintf (q, "# V=%d (%g,%g)\n", jv, x, y); printf ("%s", q), fprintf (yfile, "%s", q); } /* 572-577 */ } /* 566-578 */ } /* 564-579 */ goto Draw_Lines; /* 646 */ } /* 562-581 */ /* The two vertices in common are KV and LV. */ e = gxov(kv), f = gyov(kv), g = gxov(lv), h = gyov(lv); s = hxov(kv), t = hyov(kv), u = hxov(lv), v = hyov(lv); if (idebug) { /* 589-596 */ sprintf (q, "#\n"); printf ("%s", q), fprintf (yfile, "%s", q); sprintf (q, "# KV %d = G(%g,%g), H(%g,%g)\n", kv, e, f, s, t); printf ("%s", q), fprintf (yfile, "%s", q); sprintf (q, "# LV %d = G(%g,%g), H(%g,%g)\n", lv, g, h, u, v); printf ("%s", q), fprintf (yfile, "%s", q); } /* 589-596 */ /* Do the distances match? */ w = ((e-g)*(e-g))+((f-h)*(f-h)); z = ((s-u)*(s-u))+((t-v)*(t-v)); if ((iclose (w, z) == FALSE) || (iclose (w, 0))) { /* 602-606 */ sprintf (q, "#! Dist (%g,%g)-(%g,%g) != (%g,%g)-(%g,%g)!\n", e, f, g, h, s, t, u, v), printf ("%s", q), fprintf (yfile, "%s", q); } /* 602-606 */ c = g-e, d = h-f; r = sqrt ((c*c)+(d*d)); c /= r, d /= r; w = u-s, z = v-t; r = sqrt ((w*w)+(z*z)); w /= r, z /= r; a = (c*w)+(d*z); b = -(d*w)+(c*z); /* Figure out new vertices in plane JP. */ for (jv = 1; jv <= jnov; jv++) { /* 622-644 */ if (nopv(jp,jv) == FALSE) continue; x = gxov(jv)-e, y = gyov(jv)-f; v = (a*x)-(b*y); w = (b*x)+(a*y); hxov(jv) = v+s, hyov(jv) = w+t; if ((jv != kv) && (jv != lv)) for (mp = 1; mp <= jnop; mp++) mopv(mp,jv) = FALSE; mopv(jp,jv) = TRUE; if (idebug) { /* 637-643 */ sprintf (q, "#\n"); printf ("%s", q), fprintf (yfile, "%s", q); sprintf (q, "# V=%d (%g,%g)-->(%g,%g)\n", jv, x+e, y+f, v+s, w+t); printf ("%s", q), fprintf (yfile, "%s", q); } /* 637-643 */ } /* 622-644 */ Draw_Lines: /* 646 */ /* Any two points in plane JP that share some */ /* other plane in common get connected with a line. */ /* (If IDEBUG is on, then spread things out a bit.) */ fprintf (yfile, "#\n"); for (jv = 1; jv <= jnov; jv++) { /* 654-686 */ if (nopv(jp,jv) == FALSE) continue; u = (hxov(jv)*100.0)+(idebug*jp); v = (hyov(jv)*100.0)+(idebug*jp); for (mv = jv+1; mv <= jnov; mv++) { /* 660-685 */ if (nopv(jp,mv) == FALSE) continue; for (mp = 1; mp <= jnop; mp++) if ((mp != jp) && (nopv(mp,jv)) && (nopv(mp,mv))) break; if (mp <= jnop) { /* 668-684 */ x = (hxov(mv)*100.0)+(idebug*jp); y = (hyov(mv)*100.0)+(idebug*jp); sprintf (q, "l %d %g %g %g %g\n", jp, u, v, x, y); printf ("%s", q), fprintf (yfile, "%s", q); jnol++; if (jnol > NUML) { /* 677-680 */ printf ("Too many lines (%d)\n", jnol); exit (9); } /* 677-680 */ npol(jnol) = jp; fxol(jnol) = u, fyol(jnol) = v; gxol(jnol) = x, gyol(jnol) = y; } /* 668-684 */ } /* 660-685 */ } /* 654-686 */ /* Put letter in the center of the face. */ /* (It might be too large but put it there anyway.) */ x = y = 0.0, z = EPS; for (jv = 1; jv <= jnov; jv++) if (nopv(jp,jv)) x += hxov(jv), y += hyov(jv), z += 0.01; x = (x/z)-35.0, y = (y/z)-45.0; fprintf (yfile, "t %d %g %g 0.5 0 0 1\n", jp, x, y); fprintf (yfile, "s 3 0 0 #%c\n", qnop(jp)); fprintf (yfile, "#\n"); } /* 468-701 */ } /* 464-702 */ fclose (yfile); draw_picture (); return (0); } /* 293-709 */ /* ***** DRAW_PICTURE ***** */ void draw_picture (void) { /* 713-818 */ /* On a personal computer we can draw */ /* all the lines using the EGA mode. */ #if PC #define MAXX 640 #define MAXY 350 #define BIG_NUM 1.0e30 int ix1, iy1, ix2, iy2, icolor; int ixinc, iyinc, isflag, itemp; int idx, idy, ia, ib, itv; int maxx, maxy, maxc, jl; double x_max, x_min, y_max, y_min, f; static union REGS r; /* Are there any lines to draw? */ if (jnol == 0) { /* 730-733 */ printf ("No lines to draw; set is probably infeasible.\n"); return; } /* 730-733 */ /* Find minimum and maximum X and Y. */ x_min = BIG_NUM, x_max = -BIG_NUM; y_min = BIG_NUM, y_max = -BIG_NUM; for (jl = 1; jl <= jnol; jl++) { /* 739-744 */ x_min = min (x_min, fxol(jl)), x_max = max (x_max, fxol(jl)); y_min = min (y_min, fyol(jl)), y_max = max (y_max, fyol(jl)); x_min = min (x_min, gxol(jl)), x_max = max (x_max, gxol(jl)); y_min = min (y_min, gyol(jl)), y_max = max (y_max, gyol(jl)); } /* 739-744 */ /* Add a bit in case MIN=MAX. */ x_max += 0.001, y_max += 0.001; /* Adjust aspect ratio of screen. */ f = (y_max-y_min)*1.3/(x_max-x_min); if (f > 1.0) x_max = ((x_max-x_min)*f)+x_min; else y_max = ((y_max-y_min)/f)+y_min; /* And flip Y-axis for EGA display. */ f = y_max, y_max = y_min, y_min = f; /* Set up EGA mode. */ r.x.ax = 16, r.x.bx = 0, int86 (16, &r, &r); maxx = 639, maxy = 349, maxc = 15; for (jl = 1; jl <= jnol; jl++) { /* 759-813 */ /* First point. */ ix1 = (fxol(jl)-x_min)*maxx/(x_max-x_min); iy1 = (fyol(jl)-y_min)*maxy/(y_max-y_min); /* Second point. */ ix2 = (gxol(jl)-x_min)*maxx/(x_max-x_min); iy2 = (gyol(jl)-y_min)*maxy/(y_max-y_min); /* Color. (High numbers are brighter on EGA.) */ icolor = maxc-(npol(jl)%(maxc/2)); /* Delta-X and direction increment. */ idx = ix2-ix1, ixinc = 1; if (idx < 0) idx = -idx, ixinc = -1; /* Delta-Y and direction increment. */ idy = iy2-iy1, iyinc = 1; if (idy < 0) idy = -idy, iyinc = -1; /*IDX must be greater than or equal to IDY. */ isflag = FALSE; if (idy > idx) isflag = TRUE, itemp = idx, idx = idy, idy = itemp; /* Factor to add after each horizontal move. */ ia = idy*2; /* Test variable for the move. */ itv = ia-idx; /* Factor to add after each vertical move. */ ib = (idy-idx)*2; idx++; while (idx > 0) { /* 793-811 */ /* Put a dot on the screen. */ r.x.ax = 3072+icolor, r.x.bx = 0, r.x.cx = ix1, r.x.dx = iy1; int86 (16, &r, &r); if (itv < 0) { /* 799-801 */ itv = itv+ia; } /* 799-801 */ else { /* 803-807 */ itv = itv+ib; if (isflag) ix1 = ix1+ixinc; else iy1 = iy1+iyinc; } /* 803-807 */ if (isflag) iy1 = iy1+iyinc; else ix1 = ix1+ixinc; idx--; } /* 793-811 */ if (idebug) getch (); } /* 759-813 */ /* Go back to text mode after user types character. */ getch (), r.x.ax = 3, int86 (16, &r, &r); #endif return; } /* 713-818 */ /* ***** ICLOSE ***** */ int iclose (double a, double b) { /* 822-827 */ /* One if two values are within */ /* EPSilon and zero otherwise. */ if (((a-b < EPS) && (b-a) < EPS)) return (TRUE); else return (FALSE); } /* 822-827 */ /* ***** INVERT ***** */ void invert (void) { /* 831-850 */ det = (a*e*o)+(b*f*g)+(c*d*h)-(a*f*h)-(b*d*o)-(c*e*g); if (iclose (det, 0.0)) { /* 834-837 */ det = 0.0; return; } /* 834-837 */ aa = ((e*o)-(h*f))/det; bb = ((c*h)-(b*o))/det; cc = ((b*f)-(e*c))/det; dd = ((f*g)-(d*o))/det; ee = ((a*o)-(c*g))/det; ff = ((c*d)-(a*f))/det; gg = ((d*h)-(e*g))/det; hh = ((b*g)-(a*h))/det; oo = ((a*e)-(b*d))/det; return; } /* 831-850 */ /* ***** MAL ***** */ void *mal (int n) { /* 854-860 */ void *p; p = malloc (n); if (p == NULL) printf ("Memory allocation failure.\n"), exit (9); return (p); } /* 854-860 */ /* ***** MEMORY_ALLOCATION ***** */ void memory_allocation (void) { /* 864-903 */ int jp; /* Line arrays. */ Fxol = (double *) mal ((NUML+1)*sizeof (double)); Fyol = (double *) mal ((NUML+1)*sizeof (double)); Gxol = (double *) mal ((NUML+1)*sizeof (double)); Gyol = (double *) mal ((NUML+1)*sizeof (double)); Npol = (int *) mal ((NUML+1)*sizeof (int)); /* Plane arrays. */ Qnop = (char *) mal ((NUMP+1)*sizeof (char)); Fxop = (double *) mal ((NUMP+1)*sizeof (double)); Fyop = (double *) mal ((NUMP+1)*sizeof (double)); Fzop = (double *) mal ((NUMP+1)*sizeof (double)); Fbop = (double *) mal ((NUMP+1)*sizeof (double)); Iop = (int *) mal ((NUMP+1)*sizeof (int)); /* Vertex arrays. */ Fxov = (double *) mal ((NUMV+1)*sizeof (double)); Fyov = (double *) mal ((NUMV+1)*sizeof (double)); Fzov = (double *) mal ((NUMV+1)*sizeof (double)); Gxov = (double *) mal ((NUMV+1)*sizeof (double)); Gyov = (double *) mal ((NUMV+1)*sizeof (double)); Hxov = (double *) mal ((NUMV+1)*sizeof (double)); Hyov = (double *) mal ((NUMV+1)*sizeof (double)); /* Plane by vertex arrays. */ Nopv = (int **) mal ((NUMP+1)*sizeof (int *)); Mopv = (int **) mal ((NUMP+1)*sizeof (int *)); for (jp = 1; jp <= NUMP; jp++) { /* 898-901 */ Nopv[jp] = (int *) mal ((NUMV+1)*sizeof (int)); Mopv[jp] = (int *) mal ((NUMV+1)*sizeof (int)); } /* 898-901 */ return; } /* 864-903 */ /* Array bounds checking. In C there is no */ /* "automatic" range checker like the ones that */ /* usually come with FORTRAN or PASCAL compilers. */ /* ***** WATCH ***** */ void watch (void) { /* 911-914 */ /* This is here for debugger break points. */ return; } /* 911-914 */ /* ***** ZRANGE ***** */ int zrange (int m, int nmin, int nmax, char *o, int l, char *p) { /* 918-945 */ char z[100]; int j; static int kount = 0; static int lline = 0; if (l != lline) kount = 1, lline = l; else kount++; if ((m < nmin) || (m > nmax)) { /* 928-943 */ if (m < nmin) { /* 930-934 */ j = sprintf (z, "Value %s=%d is less than %d.", p, m, nmin); if (kount > 1) sprintf (z+j, " (%d)", kount); m = nmin; } /* 930-934 */ else if (m > nmax) { /* 936-940 */ j = sprintf (z, "Value %s=%d is greater than %d.", p, m, nmax); if (kount > 1) sprintf (z+j, " (%d)", kount); m = nmax; } /* 936-940 */ printf ("%s(%d): %s\n", o, l, z); zzrange (1); } /* 928-943 */ return (m); } /* 918-945 */ /* ***** ZZRANGE ***** */ void zzrange (int i) { /* 949-959 */ static j_err_count = 0; j_err_count += i; if (j_err_count < -1000) j_err_count = -1000; if (j_err_count > 5) { /* 954-957 */ printf ("Too many range err" "ors.\n"); exit (2); } /* 954-957 */ return; } /* 949-959 */ /* *** END OF PLANES.C *** */