/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "da3d.p" */ #include /* da3d: diana da file to 3d graphics Tom Schneider NCI/FCRDC Bldg 469. Room 144 P.O. Box B Frederick, MD 21702-1201 (301) 846-5581 (-5532 for messages) network address: toms@ncifcrf.gov National Cancer Institute Laboratory of Mathematical Biology */ /* end of program */ /* begin module version */ #define version 1.23 /* of da3d.p 1995 Nov 13 origin 1991 December 11 */ /* end module version */ /* begin module describe.da3d */ /* name da3d: diana da file to 3d graphics synopsis da3d(da: in, da3dp: in, scene: out, output: out) files da: output of the diana program; position to poistion correlations da3dp: parameter file to control scene. horizontal: shift of graph horizontally (in cm) vertical: shift of graph vertically (in cm) xlocation: location of viewer in bases ylocation: location of viewer in bases zlocation: location of viewer in bits magnify: magnification factor for whole scene, 1 = no change. xmagnify: magnification factor for x axis only. ymagnify: magnification factor for y axis only. zmagnify: magnification factor for z axis only. datacolumn: (integer) column of da which to use for the graph. dodiagonal: if the first character is 't', the diagonal values will be included. scene: 3D scene of the da data according to da3dp. Result is in PostScript output: messages to the user description Show the position to position correlation data in three dimensions. examples documentation see also diana.p author Thomas Dana Schneider bugs history [1995 Nov 13] Put lines to previous data trace so that it looks like a set of squares rather than just lines. technical notes */ /* end module describe.da3d */ /* begin module da3d.const */ #define mingrid (-100) /* minimum of side of the data grid */ #define maxgrid 100 /* maxmum of side of the data grid */ /* end module da3d.const */ /* begin module interact.const */ #define maxstring 150 /* the maximum string */ /* end module interact.const version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.const */ #define pi 3.14159265354 /* circumference divided by diameter of circle */ #define picfield 8 /* width of numbers printed to the file */ #define picwidth 5 /* number of decimal places for numbers */ #define charwidth 0.0625 /* the width of characters in inches (ie, inches/char) this allows centering of strings. */ /* note: for the Times-Roman font, 0.0625 is a good value. for the Courier-Bold font, 0.08 is a good value. */ #define dotfactor 0.00625 /* the size of dots */ #define defscale 81 /* default scale factor. coordinate units per inch */ /* end module pic.const version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.3d.type */ /* these types are used by the three dimensional graphics routines */ typedef double threevector[3]; /* a point in 3 space */ typedef double tbtarray[3][3]; /* a three by three array */ typedef struct screen { /* define a screen for viewing a 3d object */ threevector a; /* center of screen */ threevector b; /* screen x coordinate direction */ threevector c; /* screen y coordinate direction */ threevector v; /* the position of the viewer */ threevector g; /* gaze: viewing direction */ double smag; /* the magnification factor for the screen */ double range; /* 1/smag; the half width of the screen */ } screen; /* end module pic.3d.type version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module interact.type */ typedef struct string { /* a string of characters */ Char letters[maxstring]; /* the letters in the string */ long length; /* the number of characters in the string */ long current; /* the letter we are working on */ } string; /* end module interact.type version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module da3d.type */ /* datagrid = array[mingrid,maxgrid..mingrid..maxgrid] of real;*/ typedef double datagrid[maxgrid - mingrid + 1][maxgrid - mingrid + 1]; /* end module da3d.type */ Static _TEXT da; /* input file */ Static _TEXT da3dp; /* parameters */ /* PostScript output */ Static _TEXT scene; /* a file used by this program */ /* begin module pic.var */ Static boolean inpicture; /* true if we are drawing the picture, ie, startpic has been called */ Static double picxglobal, picyglobal; /* absolute location in the graph */ Static double pictolerance; /* 10 raised to the picwidth, to detect values close to zero */ Static double scale; /* scale factor. graphic coordinate units per inch */ /* NONSTANDARD for efficient use of postscript, keep track of whether there is a current path */ Static boolean inpath; /* NONSTANDARD keep track of number of segments drawn so that they can be stroked. This (probably) solves the problem of the Apple printer dying because it can't handle the data. */ Static long segments; Static double xsideold, ysideold; /* current size of a rectangle. see rectsize */ Static jmp_buf _JL1; /* end module pic.var version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module halt */ Static Void halt() { /* stop the program. the procedure performs a goto to the end of the program. you must have a label: label 1; declared, and also the end of the program must have this label: 1: end. examples are in the module libraries. this is the only goto in the delila system. */ printf(" program halt.\n"); longjmp(_JL1, 1); } /* end module halt version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module interact.clearstring */ Static Void clearstring(ribbon) string *ribbon; { /* empty the string */ long index; /* to the ribbon */ for (index = 0; index < maxstring; index++) ribbon->letters[index] = ' '; ribbon->length = 0; ribbon->current = 0; } /* clearstring */ /* end module interact.clearstring version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.functions */ /* ********************************************************************** */ /* begin module pic.await */ Static Void await() { /* Wait for user to type a carriage return. the routine assumes that there is a global file called input. */ /* the infinite way: writeln(output); writeln(output,'*********************************'); writeln(output,'* Use control-c to kill program *'); writeln(output,'*********************************'); while true do begin end;*/ } /* end module pic.await version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.startpic */ Static Void startpic(afile, setscale, x, y, thefont) _TEXT *afile; double setscale, x, y; Char thefont; { /* open the graphics field, with the given scale, and at (x,y) in that scale. scale is in device coordinates per inch. The font is chosen with thefont; t = Times-Roman, c = Courier-Bold */ /* start pic output to file afile, set the globals */ /* NONSTANDARD */ /* this is the actual "world" coordinates used: */ /* xmin, xmax, ymin, ymax */ /* ns; if (setwindow(-5.0/scale, +5.0/scale, -5.0/scale, +5.0/scale)*/ fprintf(afile->f, "gsave\n"); /* save the current graphics state */ fprintf(afile->f, "initgraphics\n"); /* make sure the printer is ready to print, without this, sometimes an Apple laserwriter will print the graph upside down, tiny and backwards! */ scale = setscale; /* set the global scale */ fprintf(afile->f, "clear erasepage\n"); /* clean residue from before */ switch (thefont) { case 'c': fprintf(afile->f, "/Courier-Bold findfont\n"); /* locate the font */ fprintf(afile->f, "%d scalefont\n", 10); /* set the font size in points*/ break; case 't': fprintf(afile->f, "/Times-Roman findfont\n"); /* locate the font */ fprintf(afile->f, "%d scalefont\n", 12); /* set the font size in points*/ break; } fprintf(afile->f, "setfont\n"); /* put the font into the current font */ /* If the following statement is done then it will work on the sun, but will kill the applewriter!!!! Sun's non-standard PostScript extension, setlinewidth has default 0, as stated in the Read This First and the NeWS Manual. This draws very quickly with 1 bit wide lines. If '1 setlinequality' is not done, then one cannot set the width of lines. So to use PostScript on the screen, I must first do '1 setlinequality'. However, if I send this code to the Applewriter, it kills PostScript on the Applewriter and I get no output whatsoever! (It took me several hours to figure this out, since once PostScript is killed on the Applewriter, the NEXT output is also smashed and I had to figure that out also...) So a standard PostScript program will not work correctly with the default. "Correcting" the PostScript program so that it works on the Sun means that it BOMBS on the Applewriter. The default for setlinequality should be '1 setlinequality' so that the same PostScript code can be used both on the Sun and with other devices. If you want speed, use the nonstandard form. An alternative is to redefine setlinequality so that '0 setlinequality' does give correct results with standard PostScript. Please review this, Randy. I think that Sun should fix it. writeln(afile,'1 setlinequality'); makes lines at least 1 bit wide */ /* set the scale to inches writeln(afile, scale:picfield:picwidth,' ', scale:picfield:picwidth,' scale'); */ /* define some things in postscript */ /* doline allows less stuff to be put in the output file. it takes two numbers off the stack, copies them, draws a line to them as coordinates. */ /* replaced by 'currentpoint translate' writeln(afile,'/doline { 2 copy lineto } def'); */ /* define a function that makes inches out of a number */ /* do this all internally here, it's faster writeln(afile,'/i { ',scale:picfield:picwidth,' mul} def'); */ /* move to the start point on the page */ fprintf(afile->f, "%*.*f %*.*f translate\n\n", picfield, picwidth, x * scale, picfield, picwidth, y * scale); fprintf(afile->f, "%% Define functions so the text produced is smaller\n"); fprintf(afile->f, "/a {stroke newpath 0 0} def %% special for arc\n"); fprintf(afile->f, "/c {stroke 0 0 moveto} def %% current point\n"); fprintf(afile->f, "/f {findfont 10 scalefont setfont} def\n"); fprintf(afile->f, " %% to set fonts simply use the f function. Example:\n"); fprintf(afile->f, " %%/Symbol f (\\142) /Courier-Bold f (-galactosidase\n"); fprintf(afile->f, "/l {lineto} def\n"); fprintf(afile->f, "/m {moveto} def\n"); fprintf(afile->f, "/n {stroke newpath 0 0 moveto} def\n"); /* new segment */ fprintf(afile->f, "/rl {rlineto} def\n"); fprintf(afile->f, "/rm {rmoveto} def\n"); fprintf(afile->f, "/s {newpath 0 0 moveto} def %% Start path \n"); fprintf(afile->f, "/t {currentpoint translate} def %% translate \n"); fprintf(afile->f, "/x {show} def %% show teXt \n\n"); /* start out the pathway */ inpath = false; /* start the number of segments written: */ segments = 0; /* now for the normal pic stuff: */ inpicture = true; picxglobal = 0.0; picyglobal = 0.0; pictolerance = (long)(exp(picwidth * log(10.0)) + 0.5); /*;writeln(output,'pictolerance = ',pictolerance:picfield:picwidth);*/ } /* end module pic.startpic version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.stoppic */ Static Void stoppic(afile) _TEXT *afile; { /* stop pic output to file afile */ /* NONSTANDARD */ if (inpath) { fprintf(afile->f, "stroke\n"); inpath = false; } fprintf(afile->f, "showpage\n"); fprintf(afile->f, "grestore\n"); /* restore the current graphics state to what it was before the startpic */ await(); inpicture = false; } #define buffer 10 Local Void checkseg(afile) _TEXT *afile; { /* NONSTANDARD checks how many segments have been written, if more than 'buffer', stroke them to the postscript page */ if (segments >= buffer) { fprintf(afile->f, "n\n"); segments = 0; } /* New segment: writeln(afile,'stroke newpath 0 0 moveto'); */ else segments++; } #undef buffer /* end module pic.stoppic version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.drawr */ Static Void drawr(afile, dx, dy, visibility, spacing) _TEXT *afile; double dx, dy; Char visibility; double spacing; { /* make a line to file afile by relative draw of dx,dy with visibility i invisible - dashed . dotted l line with the dashes or dots separated by the spacing given (this has no effect with invisible and line). */ /* NONSTANDARD */ double ddx, ddy; /* changes in dx and dy for dots and dashes */ double dr; /* the hypotenuse, the distance actually drawn */ boolean on; /* draw linesegment if true */ double y; /* the variable for tracking dots and dashes */ long r; /* number of times to cycle for dots and dashes */ double ss; /* precalculated value to make things a bit faster */ double theta; /* angle of the line */ long FORLIM; if (!inpath) { fprintf(afile->f, "s\n"); inpath = true; } /* starts from current coordinates */ /* Start path: writeln(afile,'newpath 0 0 moveto'); */ else checkseg(afile); /* checks if not (visibility in ['l','i','.','-']) then writeln(afile,'%YELLLLLL!!!',visibility,'!'); writeln(afile,'% ',visibility,' line');*/ /* put these on the stack, they will always be used */ fprintf(afile->f, "%*.*f %*.*f", picwidth, picfield, dx * scale, picwidth, picfield, dy * scale); switch (visibility) { case 'l': case 'i': switch (visibility) { case 'i': fprintf(afile->f, " m"); break; case 'l': fprintf(afile->f, " l"); break; } break; case '.': case '-': /* make up our own dots and dashes */ putc('\n', afile->f); /* move away from the (dx,dy) on the stack */ if (spacing <= 0.0) { printf("drawr: spacing zero with . or - line\n"); halt(); } if (dx == 0.0) { ddx = 0.0; /* avoid division by zero */ ddy = scale * spacing; if (dy < 0) ddy = -ddy; /* this makes sure that we draw lines straight down if that was the request */ } else { /* find out the angle of the slope, intentionally lose the sign */ theta = atan(fabs(dy / dx)); ddx = scale * spacing * cos(theta); ddy = scale * spacing * sin(theta); /* return the sign to the little buggers */ if (dx < 0) ddx = -ddx; if (dy < 0) ddy = -ddy; } y = 0.0; switch (visibility) { case '.': ss = scale * dotfactor; break; case '-': on = true; break; } dr = sqrt(dx * dx + dy * dy); FORLIM = (long)floor(dr / spacing + 0.5); for (r = 1; r <= FORLIM; r++) { switch (visibility) { case '-': fprintf(afile->f, "%*.*f %*.*f", picwidth, picfield, ddx, picwidth, picfield, ddy); if (on) fprintf(afile->f, " rl\n"); else fprintf(afile->f, " rm\n"); on = !on; break; case '.': fprintf(afile->f, "%*.*f 0 rl", picwidth, picfield, ss); fprintf(afile->f, " %*.*f 0 rl", picwidth, picfield, -ss); fprintf(afile->f, " %*.*f %*.*f", picwidth, picfield, ddx, picwidth, picfield, ddy); fprintf(afile->f, " rm\n"); break; /* put out a dot like in dotr */ } } /* let's make really sure we got there!! */ fprintf(afile->f, " m\n"); /* pulled from the stack */ break; } /* an elegant way to make postscript keep a global record is to translate the coordinates! */ /* writeln(afile,' currentpoint translate'); */ fprintf(afile->f, " t\n"); picxglobal += dx; picyglobal += dy; } /* end module pic.drawr version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.mover */ Static Void mover(afile, dx, dy) _TEXT *afile; double dx, dy; { /* move relative the amount (dx, dy). */ drawr(afile, dx, dy, 'i', 0.0); } /* end module pic.mover version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.liner */ Static Void liner(afile, dx, dy) _TEXT *afile; double dx, dy; { /* draw a line the relative amount (dx, dy). */ drawr(afile, dx, dy, 'l', 0.0); } /* end module pic.liner version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.drawa */ Static Void drawa(afile, x, y, visibility, spacing) _TEXT *afile; double x, y; Char visibility; double spacing; { /* make a line to file afile to absolute coordinate x,y with visibility i invisible - dashed . dotted l line with the dashes or dots separated by the spacing given (this has no effect with invisible and line). */ double dx, dy; /* differences between current and desired locations */ dx = x - picxglobal; dy = y - picyglobal; drawr(afile, dx, dy, visibility, spacing); } /* end module pic.drawa version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.movea */ Static Void movea(afile, x, y) _TEXT *afile; double x, y; { /* move to absolute x and y */ drawa(afile, x, y, 'i', 0.0); } /* end module pic.movea version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.linea */ Static Void linea(afile, x, y) _TEXT *afile; double x, y; { /* draw a line from current position to absolute x and y */ drawa(afile, x, y, 'l', 0.0); } /* end module pic.linea version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.graphstring */ Static Void graphstring(tofile, s, centered) _TEXT *tofile; string *s; boolean centered; { /* graph the string s. If it is recognized as a quoted string (surrounded by double quotes), graph it without the quotes and center it. Always center if centered is true. Otherwise simply graph it. if not in picture, just write it to output */ /* NONSTANDARD */ long i; /* index to s, and temporary storage */ double mv; /* holds amount to move, in plotting coordinates */ boolean quoted; /* true if the string is quoted */ long FORLIM; if (!(inpicture && s->length > 0)) return; if (s->length > 2) { if (s->letters[0] == '"' && s->letters[s->length - 1] == '"') quoted = true; else quoted = false; } else quoted = false; /* override so quoted strings are always centered */ if (quoted) centered = true; if (centered) { /* generate the calls to center the string. Note: this must not be done for the pic program, which already centers */ if (quoted) i = s->length - 2; else i = s->length; mv = i * charwidth / 2.0; mover(tofile, -mv, 0.0); } /* do the non-standard postscript: */ /* do postscript to complete pervious path */ /* set current point: writeln(tofile,'stroke 0 0 moveto'); */ fprintf(tofile->f, "c\n"); putc('(', tofile->f); /* begin postscript literal */ if (quoted) { FORLIM = s->length - 2; for (i = 1; i <= FORLIM; i++) putc(s->letters[i], tofile->f); } else { FORLIM = s->length; for (i = 0; i < FORLIM; i++) putc(s->letters[i], tofile->f); } putc(')', tofile->f); /* end postscript literal */ fprintf(tofile->f, " x\n"); /* show the literal */ inpath = false; /* force new path from here */ if (centered) { /* restore to previous location */ mover(tofile, mv, 0.0); /* no output if not in picture writestring(tofile,s); writeln(tofile) */ } } /* end module pic.graphstring version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.stringinteger */ Static Void stringinteger(number, name, width, leadingzeros) long number; string *name; long width; boolean leadingzeros; { /* make the string from the number, start putting characters in after the current length point. use width characters. if leadingzeros is true, trail zeros before the number. */ long bigdigit; /* the location of the biggest digit */ long dig; /* number of digits in the number */ long place; /* place to write the next digit of the number */ long sign; /* the sign of the number */ if (number < 0) { sign = -1; name->length++; /* provide room for the sign!! */ number = -number; if (leadingzeros) printf( "WARNING: stringinteger: the sign of a negative number with leading zeros is lost\n"); } else sign = 1; /* log 10 of the number plus 1 is the number of digits in the number. On this sun computer ln(1000)/ln(10) is 2.9999, which when truncated gives 2, rather than the desired 3. To avoid this kind of problem, 0.1 is added. */ if (number > 9) dig = (long)(log(number + 0.1) / log(10.0)) + 1; else dig = 1; if (dig > width) { printf("stringinteger: number width too small\n"); printf("%ld digit number (%ld)\n", dig, number); printf("does not fit in %ld characters\n", width); halt(); } if (leadingzeros) bigdigit = name->length + 1; /* no sign if leading zeros */ else { bigdigit = name->length + width - dig + 1; if (bigdigit <= name->length && sign < 0) { printf("stringinteger: no room for sign\n"); halt(); } } if (sign < 0) name->letters[bigdigit-2] = '-'; for (place = name->length + width - 1; place >= bigdigit - 1; place--) { /* p2c: da3d.p, line 632: * Note: Using % for possibly-negative arguments [317] */ switch (number % 10) { case 0: name->letters[place] = '0'; break; case 1: name->letters[place] = '1'; break; case 2: name->letters[place] = '2'; break; case 3: name->letters[place] = '3'; break; case 4: name->letters[place] = '4'; break; case 5: name->letters[place] = '5'; break; case 6: name->letters[place] = '6'; break; case 7: name->letters[place] = '7'; break; case 8: name->letters[place] = '8'; break; case 9: name->letters[place] = '9'; break; } number /= 10; } name->length += width; } /* end module pic.stringinteger version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.stringreal */ Static Void stringreal(number, name, width, decimal) double number; string *name; long width, decimal; { /* make the string from the real number, start putting characters in at the start point. use width characters and decimal characters after the decimal place */ /* note that the rounding operation to get the digits below zero must be done first. then the digits above zero can be lopped off. this makes 99.99 come out correctly to 100.0 (to 1 decimal place) otherwise, 99.99 -> 0.99 -> 1.0 (rounded) -> 10 (print with 1 decimal place), and stringinteger won't be happy about that. */ long abovezero; /* the number shifted above the decimal place, to 'decimal' positions (and rounded) */ long shift; /* power of ten used to shift a number around relative to the decimal point */ long sign; /* the sign of the number */ long thedecimal; /* integer version of the decimal part of the number */ long theupper; /* integer version of the upper part of the number */ if (number < 0) sign = -1; else sign = 1; number = fabs(number); /* make positive */ /* the amount to shift the number above zero */ shift = (long)floor(exp(decimal * log(10.0)) + 0.5); /* amount to move above zero */ abovezero = (long)floor(number * shift + 0.5); /* move above zero, round off */ theupper = (long)((double)abovezero / shift); thedecimal = abovezero - shift * theupper; /* create the actual real number */ /* before decimal point */ stringinteger(sign * theupper, name, width - decimal - 1, false); /* put in the decimal point */ name->length++; name->letters[name->length - 1] = '.'; stringinteger(thedecimal, name, decimal, true); /* after decimal point */ } /* end module pic.stringreal version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.picnumber */ Static Void picnumber(afile, dx, dy, number, width, decimal, centered) _TEXT *afile; double dx, dy, number; long width, decimal; boolean centered; { /* Supply graphic commands for a 'number' whose center is at the relative point (dx, dy) from the current point, 'width' characters wide and 'decimal' characters beyond the decimal point. If the width is zero, no number is produced. procedure stringnumber(number: integer; start: integer; var name: string); the location after the call is the same as before the call. The string is optionally centered */ string name; /* the string to pack the number into for shipping out */ if (width <= 0) return; mover(afile, dx, dy); clearstring(&name); if (decimal > 0) stringreal(number, &name, width, decimal); else stringinteger((long)floor(number + 0.5), &name, width, false); graphstring(afile, &name, centered); mover(afile, -dx, -dy); } /* end module pic.picnumber version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.xtic */ Static Void xtic(afile, length, dx, dy, number, width, decimal) _TEXT *afile; double length, dx, dy, number; long width, decimal; { /* produce a tic mark for the x axis of "length" long. Supply a number whose center is at the relative point (dx, dy) from the end to the tick, 'width' characters wide and 'decimal' characters beyond the decimal point. If the width is zero, no number is produced. the location after the call is the same as before the call. */ liner(afile, 0.0, -length); picnumber(afile, dx, dy, number, width, decimal, true); mover(afile, 0.0, length); } /* end module pic.xtic version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.ytic */ Static Void ytic(afile, length, dx, dy, number, width, decimal) _TEXT *afile; double length, dx, dy, number; long width, decimal; { /* produce a tic mark for the y axis of "length" long. Supply a number whose center is at the relative point (dx, dy) from the end to the tick, 'width' characters wide and 'decimal' characters beyond the decimal point. If the width is zero, no number is produced. the location after the call is the same as before the call. */ liner(afile, -length, 0.0); picnumber(afile, dx, dy, number, width, decimal, true); mover(afile, length, 0.0); } /* end module pic.ytic version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.xaxis */ Static Void xaxis(afile, axlength, fromtic, interval, totic, length, dx, dy, width, decimal) _TEXT *afile; double axlength, fromtic, interval, totic, length, dx, dy; long width, decimal; { /* draw an x axis starting from the current position. the length of the xaxis is axlength. the axis is labeled with numbers starting with fromtic at intervals given up to totic. the remaining variables describe the form of the tic marks as in xtic. If the width is zero, no number is produced. the location after the call is the same as before the call. */ double jump; /* the space to move on the graph between tic marks */ double jumpdistance = 0.0; /* the total jumps made. this may not be a simple function of the input variables since they may not work out to an exact number of jumps */ double tic; /* the numerical value of the tic label */ double totici2; /* an extra half interval guarantees the last tic is placed */ liner(afile, axlength, 0.0); mover(afile, -axlength, 0.0); if (totic == fromtic) { printf("xaxis: fromtic and totic cannot be equal\n"); halt(); } if (axlength == 0.0 || interval == 0.0) { printf("xaxis: neither axlength nor interval can be zero\n"); halt(); } jump = axlength * interval / (totic - fromtic); /* note! the extra amount of interval/2 makes sure that the last tic is placed, if the computer does not add quite right and overshoots by a tich on the value of tic. This much makes it surely work! */ tic = fromtic; if (interval > 0.0) { totici2 = totic + interval / 2.0; while (tic <= totici2) { xtic(afile, length, dx, dy, tic, width, decimal); tic += interval; if (tic > totici2) break; mover(afile, jump, 0.0); jumpdistance += jump; } } else if (interval < 0.0) { totici2 = totic - interval / 2.0; while (tic >= totici2) { xtic(afile, length, dx, dy, tic, width, decimal); tic += interval; if (tic < totici2) break; mover(afile, jump, 0.0); jumpdistance += jump; } } mover(afile, -jumpdistance, 0.0); } /* end module pic.xaxis version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.yaxis */ Static Void yaxis(afile, aylength, fromtic, interval, totic, length, dx, dy, width, decimal) _TEXT *afile; double aylength, fromtic, interval, totic, length, dx, dy; long width, decimal; { /* draw a y axis starting from the current position. the length of the yaxis is aylength. the axis is labeled with numbers starting with fromtic at intervals given up to totic. the remaining variables describe the form of the tic marks as in ytic. If the width is zero, no number is produced. the location after the call is the same as before the call. */ double half; /* half the jump interval. By adding this to the while loops, we assure that the very last tic gets done, and isn't lost due to roundoff */ double jump; /* the space to move on the graph between tic marks */ double jumpdistance = 0.0; /* the total jumps made. this may not be a simple function of the input variables since they may not work out to an exact number of jumps */ double tic; /* the numerical value of the tic label */ liner(afile, 0.0, aylength); mover(afile, 0.0, -aylength); if (totic == fromtic) { printf("yaxis: fromtic and totic cannot be equal\n"); halt(); } if (aylength == 0.0 || interval == 0.0) { printf("yaxis: neither aylength nor interval can be zero\n"); halt(); } jump = aylength * interval / (totic - fromtic); half = interval / 2.0; tic = fromtic; if (interval > 0.0) { while (tic <= totic + half) { ytic(afile, length, dx, dy, tic, width, decimal); tic += interval; if (tic > totic + half) break; mover(afile, 0.0, jump); jumpdistance += jump; } } else if (interval < 0.0) { while (tic >= totic - half) { ytic(afile, length, dx, dy, tic, width, decimal); tic += interval; if (tic < totic - half) break; mover(afile, 0.0, jump); jumpdistance += jump; } } mover(afile, 0.0, -jumpdistance); } /* end module pic.yaxis version = 2.63; (@ of dops.p 1991 November 2 */ /* ********************************************************************** */ /* end module pic.functions version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.polrec */ Static Void polrec(r, theta, x, y) double r, theta, *x, *y; { /* convert polar to rectangular coordinates, theta is in radians */ *x = r * cos(theta); *y = r * sin(theta); } /* end module pic.polrec version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.3d.package */ /* ********************************************************************** */ /* begin module pic.3d.determinant */ Static double determinant(a) double (*a)[3]; { /* compute the determinant of a */ return (a[0][0] * (a[1][1] * a[2][2] - a[2][1] * a[1][2]) + a[0] [1] * (a[2][0] * a[1][2] - a[1][0] * a[2][2]) + a[0] [2] * (a[1][0] * a[2][1] - a[2][0] * a[1][1])); } /* end module pic.3d.determinant version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.3d.d32 */ Static Void d32(o, a, b, c, v, xloc, yloc) double *o, *a, *b, *c, *v; double *xloc, *yloc; { /* convert from 3d to 2d. the players are: o: the coordinate of the object point to be converted to 2d a,b,c: define the position of the window (screen): a: center of screen b: screen x coordinate direction c: screen y coordinate direction v: the position of the viewer xloc,yloc: the resulting image vector in screen coordinates. The method of graphics is to project the object (o) toward the viewer (v) and to determine the interception of this line with the screen as defined by a,b and c. the result is expressed in the coordinate system of the screen, and so can be plotted on a 2d plotting device. When one works through the vector math, it turns out that to find the screen coordinates requires solving a set of linear equations. This is done using Cramer's rule and determinants. */ double ov, oa; /* for partial calculation */ long j; /* index to the arrays */ tbtarray d, x, y; /* define the coefficients of the equations in d,x and y */ for (j = 0; j <= 2; j++) { ov = o[j] - v[j]; d[j][0] = b[j]; d[j][1] = c[j]; d[j][2] = ov; oa = o[j] - a[j]; x[j][0] = oa; x[j][1] = c[j]; x[j][2] = ov; y[j][0] = b[j]; y[j][1] = oa; y[j][2] = ov; } /* use cramer's rule to find the solution */ *xloc = determinant(x) / determinant(d); *yloc = determinant(y) / determinant(d); } /* end module pic.3d.d32 version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.3d.view */ Static Void view(v, gaze, smag, a, b, c) double *v, *gaze; double smag; double *a, *b, *c; { /* this routine converts a viewing position (v) and a viewing direction (gaze), into the a,b,c values of a vertically oriented screen (ie, the screen is right side up). a is the center of the screen, b is the x axis, c is the y axis on the screen. This saves the user the trouble to make sure that b, c and the direction of viewing are orthogonal. one may magnify the view by making smag greater than one, or one may shrink the view by making smag less than one. if the viewing direction vector is not large enough, then the program halts. note: gaze is automatically converted to a unit vector. */ double db; /* magnitude of db */ double dgaze; /* magnitude of gaze */ long j; /* index to the arrays */ /* first check out the gaze direction */ dgaze = sqrt(gaze[0] * gaze[0] + gaze[1] * gaze[1] + gaze[2] * gaze[2]); if (smag == 0.0) { printf("screen magnitude cannot be zero\n"); halt(); } if (dgaze <= 0.001) { printf("gaze magnitude (%5.3f) is too small\n", dgaze); halt(); } /* make gaze a unit vector and set up the a vector as the viewing point plus the gaze vector */ for (j = 0; j <= 2; j++) { gaze[j] /= dgaze; a[j] = v[j] + gaze[j]; } /* the x axis of the screen, the b vector, is horizontal and orthogonal to the gaze */ b[0] = gaze[1]; b[1] = -gaze[0]; b[2] = 0.0; db = sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]); /* check for top view case and correct if so: */ if (db == 0.0) { db = 1.0; b[0] = 1.0; b[1] = 0.0; /* b[3] := 0; already from above */ } else { for (j = 0; j <= 2; j++) b[j] /= db; } /* make b a unit vector */ /* now that the gaze is a unit vector, and we have constructed the x axis in the b vector also as a unit vector, the cross product of these two will generate the y axis as a unit vector, c: */ c[0] = b[1] * gaze[2] - gaze[1] * b[2]; c[1] = gaze[0] * b[2] - b[0] * gaze[2]; c[2] = b[0] * gaze[1] - gaze[0] * b[1]; /* now normalize both b and c vectors to be of size 1/smag */ for (j = 0; j <= 2; j++) { b[j] /= smag; c[j] /= smag; } } /* end module pic.3d.view version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.3d.makescreen */ Static Void makescreen(vx, vy, vz, gx, gy, gz, smagnitude, s) double vx, vy, vz, gx, gy, gz, smagnitude; screen *s; { /* create the screen s based on the viewing location (vx,vy,vz) and the direction of gaze (gz,gy,gz). The screen size is scaled by smagnitude; doubling smagnitude should double the size of the scene. */ /* This routine makes creation of the screen very simple for the user. One need not look at the view routine. */ s->v[0] = vx; s->v[1] = vy; s->v[2] = vz; s->g[0] = gx; s->g[1] = gy; s->g[2] = gz; view(s->v, s->g, smagnitude, s->a, s->b, s->c); s->smag = smagnitude; s->range = 1 / smagnitude; } /* end module pic.3d.makescreen version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.3d.project3d */ Static Void project3d(x, y, z, s, xscreen, yscreen) double x, y, z; screen s; double *xscreen, *yscreen; { /* project the point (x,y,z) onto the screen s, to find the screen coordinates (xscreen and yscreen). */ /* This routine simplifies the projection function for the user. */ threevector o; /* for passing the values to d32 */ o[0] = x; o[1] = y; o[2] = z; d32(o, s.a, s.b, s.c, s.v, xscreen, yscreen); } /* end module pic.3d.project3d version = 2.63; (@ of dops.p 1991 November 2 */ /* ********************************************************************** */ /* end module pic.3d.package version = 2.63; (@ of dops.p 1991 November 2 */ /* ********************************************************************** */ /* ********************************************************************** */ /* ********************************************************************** */ /* begin module pic.3d.test.fun */ Static double fun(r) double r; { /* a function to plot */ return (3 / (1 + r * r / 2)); } /* end module pic.3d.test.fun version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module pic.3d.test.test3d */ Static Void test3d(afile) _TEXT *afile; { /* test three dimensional graphics */ screen s; /* the screen on which to project the 3d image */ double xscreen, yscreen; /* location on the screen corresponding to the projection of o onto the screen defined by v,a,b,c */ double oldxscreen, oldyscreen; /* the previous valuse of xscreen and yscreen */ /* definition of a spiral */ double dr; /* change in r */ double dtheta; /* change in theta */ double r = 0.0; /* radius of the current position */ double radius = 2.0; /* the radius of the spiral */ double theta = 0.0; /* angle of the current position */ double thickness = 0.1; /* spacing between spiral arms */ double steps = 15.0; /* number of steps around a circle of the spiral */ double x = 0.0, y = 0.0; double z; /* the location in three space */ makescreen(5.0, 5.0, 5.0, -1.0, -1.0, -1.0, 5.0, &s); dr = thickness / steps; dtheta = 2 * pi / steps; z = fun(r); project3d(x, y, z, s, &oldxscreen, &oldyscreen); mover(afile, oldxscreen, oldyscreen); /* premove to the startpoint of the graph */ while (r < radius) { r += dr; theta += dtheta; polrec(r, theta, &x, &y); z = fun(r); project3d(x, y, z, s, &xscreen, &yscreen); /* draw a line from where we where to the new place */ liner(afile, xscreen - oldxscreen, yscreen - oldyscreen); oldxscreen = xscreen; oldyscreen = yscreen; } } /* end module pic.3d.test.test3d version = 2.63; (@ of dops.p 1991 November 2 */ /* begin module da3d.getfromto */ Static Void getfromto(da, f, t) _TEXT *da; long *f, *t; { /* get the from (f) and to (t) range of the da file */ Char c; /* for skipping the * on the line */ if (*da->name != '\0') { if (da->f != NULL) da->f = freopen(da->name, "r", da->f); else da->f = fopen(da->name, "r"); } else rewind(da->f); if (da->f == NULL) _EscIO2(FileNotFound, da->name); RESETBUF(da->f, Char); fscanf(da->f, "%*[^\n]"); getc(da->f); /* skip id line */ fscanf(da->f, "%*[^\n]"); getc(da->f); /* skip book line */ fscanf(da->f, "%c%ld%ld%*[^\n]", &c, f, t); getc(da->f); /* skip book line */ if (c == '\n') c = ' '; printf("%ld %ld\n", *f, *t); } /* end module da3d.getfromto */ /* begin module da3d.readdata */ Static Void readdata(da, datacolumn, fromrange, torange, data) _TEXT *da; long datacolumn, *fromrange, *torange; double (*data)[maxgrid - mingrid + 1]; { /* read x,y (columns 1 and 2) and z data from column datacolumn of file da. */ Char idcharacter; /* identification character for reading data */ Char skipc; /* for skipping characters */ long column; /* for reading through columns */ long x, y; /* location on the base plane */ double z; /* the location in three space */ /* obtain the from-to range */ getfromto(da, fromrange, torange); if (*torange - *fromrange > maxgrid - mingrid) { printf( "range (from %ld to %ld) is bigger than allowed by constant maxgrid (%ld)\n", *fromrange, *torange, (long)maxgrid); halt(); } /* clear grid */ for (x = mingrid; x <= maxgrid; x++) { for (y = mingrid; y <= maxgrid; y++) data[x - mingrid][y - mingrid] = -LONG_MAX; } while (!BUFEOF(da->f)) { idcharacter = getc(da->f); if (idcharacter == '\n') idcharacter = ' '; if (idcharacter != 'i') { fscanf(da->f, "%*[^\n]"); getc(da->f); continue; } skipc = getc(da->f); skipc = getc(da->f); skipc = getc(da->f); /* skip characters */ if (skipc == '\n') skipc = ' '; if (skipc == '\n') skipc = ' '; if (skipc == '\n') skipc = ' '; fscanf(da->f, "%ld%ld", &x, &y); /* skip numbers until the z is read: */ for (column = 5; column <= datacolumn; column++) fscanf(da->f, "%lg", &z); fscanf(da->f, "%*[^\n]"); getc(da->f); if (x < mingrid || x > maxgrid || y < mingrid || y > maxgrid) { printf("point (%ld,%ld) is out of range%ld to %ld\n", x, y, (long)mingrid, (long)maxgrid); halt(); } /*zzz*/ data[x - mingrid][y - mingrid] = z; } } /* Local variables for grid: */ struct LOC_grid { _TEXT *afile; datagrid data; boolean dodiagonal; double xmagnify, ymagnify, zmagnify; screen s; /* the screen on which to project the 3d image */ double xscreen, yscreen; /* location on the screen corresponding to the projection of o onto the screen defined by v,a,b,c */ double oldxscreen, oldyscreen; /* the previous valuse of xscreen and yscreen */ long x, y; /* location on the base plane */ double z; /* the location in three space */ double xold, yold; /* the old location */ boolean changed; /* if true, the data has changed suddenly; use a move to the current coordinate rather than drawing a line */ double xhold, yhold; /* screen coordinates for temporary markings */ } ; Local Void drawit(LINK) struct LOC_grid *LINK; { /* draw to the point x,y,z */ LINK->z = LINK->data[LINK->x - mingrid][LINK->y - mingrid]; if (LINK->z <= -LONG_MAX) return; if (!(LINK->dodiagonal && LINK->x == LINK->y || LINK->x != LINK->y)) return; fprintf(LINK->afile->f, "%% (%3ld, %3ld, %10.5f)\n", LINK->x, LINK->y, LINK->z); /* any jump in coordinates restarts the drawing line */ if (fabs(LINK->x - LINK->xold) > 1 || fabs(LINK->y - LINK->yold) > 1) LINK->changed = true; else LINK->changed = false; if (LINK->changed) { /* label axis */ project3d(LINK->x * LINK->xmagnify, LINK->y * LINK->ymagnify, 0.0, LINK->s, &LINK->xhold, &LINK->yhold); movea(LINK->afile, LINK->xhold, LINK->yhold); project3d(LINK->x * LINK->xmagnify, LINK->y * LINK->ymagnify, -1.0, LINK->s, &LINK->xhold, &LINK->yhold); linea(LINK->afile, LINK->xhold, LINK->yhold); fprintf(LINK->afile->f, "gsave s (%ld) x grestore\n", LINK->x); } project3d(LINK->x * LINK->xmagnify, LINK->y * LINK->ymagnify, LINK->z * LINK->zmagnify, LINK->s, &LINK->xscreen, &LINK->yscreen); if (LINK->changed) movea(LINK->afile, LINK->xscreen, LINK->yscreen); else liner(LINK->afile, LINK->xscreen - LINK->oldxscreen, LINK->yscreen - LINK->oldyscreen); LINK->oldxscreen = LINK->xscreen; LINK->oldyscreen = LINK->yscreen; LINK->xold = LINK->x; LINK->yold = LINK->y; LINK->changed = false; } /* end module da3d.readdata */ /* begin module da3d.grid */ Static Void grid(afile_, data_, dodiagonal_, fromrange, torange, xlocation, ylocation, zlocation, magnify, xmagnify_, ymagnify_, zmagnify_) _TEXT *afile_; double (*data_)[maxgrid - mingrid + 1]; boolean dodiagonal_; double fromrange, torange, xlocation, ylocation, zlocation, magnify, xmagnify_, ymagnify_, zmagnify_; { /* create three dimensional graphics from the data */ struct LOC_grid V; V.afile = afile_; memcpy(V.data, data_, sizeof(datagrid)); V.dodiagonal = dodiagonal_; V.xmagnify = xmagnify_; V.ymagnify = ymagnify_; V.zmagnify = zmagnify_; makescreen(xlocation, ylocation, zlocation, -xlocation, -ylocation, -zlocation, magnify, &V.s); /* outline the area */ project3d(fromrange * V.xmagnify, fromrange * V.ymagnify, 0.0, V.s, &V.xscreen, &V.yscreen); movea(V.afile, V.xscreen, V.yscreen); project3d(fromrange * V.xmagnify, torange * V.ymagnify, 0.0, V.s, &V.xscreen, &V.yscreen); linea(V.afile, V.xscreen, V.yscreen); project3d(fromrange * V.xmagnify, torange * V.ymagnify, 2.0, V.s, &V.xhold, &V.yhold); linea(V.afile, V.xhold, V.yhold); fprintf(V.afile->f, "gsave s (%4.2f bits) x grestore\n", 2 / (V.zmagnify * magnify)); linea(V.afile, V.xscreen, V.yscreen); project3d(torange * V.xmagnify, torange * V.ymagnify, 0.0, V.s, &V.xscreen, &V.yscreen); linea(V.afile, V.xscreen, V.yscreen); project3d(fromrange * V.xmagnify, fromrange * V.ymagnify, 0.0, V.s, &V.xscreen, &V.yscreen); linea(V.afile, V.xscreen, V.yscreen); V.xold = -LONG_MAX; V.yold = -LONG_MAX; for (V.x = mingrid; V.x <= maxgrid; V.x++) { for (V.y = mingrid; V.y <= maxgrid; V.y++) drawit(&V); } V.changed = true; for (V.y = mingrid; V.y <= maxgrid; V.y++) { for (V.x = mingrid; V.x <= maxgrid; V.x++) drawit(&V); } } /* Local variables for themain: */ struct LOC_themain { Char dd; /* character for reading value for setting dodiagonal */ boolean dodiagonal; /* do the diagonal in the graph */ double magnify; /* factor by which to magnify the scene */ double xmagnify; /* factor by which to magnify zscale */ double ymagnify; /* factor by which to magnify zscale */ double zmagnify; /* factor by which to magnify zscale */ double xlocation; /* x location */ double ylocation; /* y location */ double zlocation; /* z location */ long datacolumn; /* column of da to use */ double horizontal; /* zero coordinate of screen in cm */ double vertical; /* zero coordinate of screen in cm */ } ; Local Void readda3d(f, LINK) _TEXT *f; struct LOC_themain *LINK; { /* obtain parameters */ fscanf(f->f, "%lg%*[^\n]", &LINK->horizontal); getc(f->f); fscanf(f->f, "%lg%*[^\n]", &LINK->vertical); getc(f->f); fscanf(f->f, "%lg%*[^\n]", &LINK->xlocation); getc(f->f); fscanf(f->f, "%lg%*[^\n]", &LINK->ylocation); getc(f->f); fscanf(f->f, "%lg%*[^\n]", &LINK->zlocation); getc(f->f); fscanf(f->f, "%lg%*[^\n]", &LINK->magnify); getc(f->f); fscanf(f->f, "%lg%*[^\n]", &LINK->xmagnify); getc(f->f); fscanf(f->f, "%lg%*[^\n]", &LINK->ymagnify); getc(f->f); fscanf(f->f, "%lg%*[^\n]", &LINK->zmagnify); getc(f->f); fscanf(f->f, "%ld%*[^\n]", &LINK->datacolumn); getc(f->f); if (LINK->datacolumn <= 4 || LINK->datacolumn > 15) { printf("data column must be >4 and < 16\n"); halt(); } fscanf(f->f, "%c%*[^\n]", &LINK->dd); getc(f->f); if (LINK->dd == '\n') LINK->dd = ' '; LINK->dodiagonal = (LINK->dd == 't'); } Local Void writeda3d(f, LINK) _TEXT *f; struct LOC_themain *LINK; { /* obtain parameters */ fprintf(f->f, "%*.*f horizontal adjust (cm)\n", picfield, picwidth, LINK->horizontal); fprintf(f->f, "%*.*f vertical adjust (cm)\n", picfield, picwidth, LINK->vertical); fprintf(f->f, "% .*E xlocation (bases)\n", P_max(picfield - 7, 1), LINK->xlocation); fprintf(f->f, "% .*E ylocation (bases)\n", P_max(picfield - 7, 1), LINK->ylocation); fprintf(f->f, "% .*E zlocation (bases)\n", P_max(picfield - 7, 1), LINK->zlocation); fprintf(f->f, "%*.*f magnify\n", picfield, picwidth, LINK->magnify); fprintf(f->f, "%*.*f xmagnify\n", picfield, picwidth, LINK->xmagnify); fprintf(f->f, "%*.*f ymagnify\n", picfield, picwidth, LINK->ymagnify); fprintf(f->f, "%*.*f zmagnify\n", picfield, picwidth, LINK->zmagnify); fprintf(f->f, "%*ld data column\n", picfield, LINK->datacolumn); if (LINK->datacolumn <= 4 || LINK->datacolumn > 15) { printf("data column must be >4 and < 16\n"); halt(); } fprintf(f->f, "%c dodiagonal\n", LINK->dd); } /* end module da3d.grid */ /******************************************************************************/ /* begin module da3d.themain */ Static Void themain(da, da3dp, scene) _TEXT *da, *da3dp, *scene; { /* the main procedure of the program */ struct LOC_themain V; datagrid data; /* store z values in x,y grid */ long fromrange, torange; /* range of the data */ double cmperinch; /* conversion factor from cm to inches */ _TEXT TEMP; printf("da3d %4.2f\n", version); if (*scene->name != '\0') { if (scene->f != NULL) scene->f = freopen(scene->name, "w", scene->f); else scene->f = fopen(scene->name, "w"); } else { if (scene->f != NULL) rewind(scene->f); else scene->f = tmpfile(); } if (scene->f == NULL) _EscIO2(FileNotFound, scene->name); SETUPBUF(scene->f, Char); fprintf(scene->f, "%%! da3d %4.2f\n\n", version); if (*da3dp->name != '\0') { if (da3dp->f != NULL) da3dp->f = freopen(da3dp->name, "r", da3dp->f); else da3dp->f = fopen(da3dp->name, "r"); } else { rewind(da3dp->f); } if (da3dp->f == NULL) _EscIO2(FileNotFound, da3dp->name); RESETBUF(da3dp->f, Char); if (BUFEOF(da3dp->f)) { test3d(scene); /* test routine */ startpic(scene, 81.0, 3.0, 5.0, 'c'); } else { readda3d(da3dp, &V); TEMP.f = stdout; *TEMP.name = '\0'; writeda3d(&TEMP, &V); cmperinch = 2.54; startpic(scene, 72.0, V.horizontal / cmperinch, V.vertical / cmperinch, 'c'); readdata(da, V.datacolumn, &fromrange, &torange, data); grid(scene, data, V.dodiagonal, (double)fromrange, (double)torange, V.xlocation, V.ylocation, V.zlocation, V.magnify, V.xmagnify, V.ymagnify, V.zmagnify); } stoppic(scene); } /* end module da3d.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; scene.f = NULL; strcpy(scene.name, "scene"); da3dp.f = NULL; strcpy(da3dp.name, "da3dp"); da.f = NULL; strcpy(da.name, "da"); themain(&da, &da3dp, &scene); _L1: if (da.f != NULL) fclose(da.f); if (da3dp.f != NULL) fclose(da3dp.f); if (scene.f != NULL) fclose(scene.f); exit(EXIT_SUCCESS); } /* End. */