/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "rsgra.p" */ #include /* rsequence graph Dr. Thomas D. Schneider National Institutes of Health National Cancer Institute Center for Cancer Research Nanobiology Program Molecular Information Theory Group Frederick, Maryland 21702-1201 toms@ncifcrf.gov permanent email: toms@alum.mit.edu (use only if first address fails) http://www.ccrnp.ncifcrf.gov/~toms/ module libraries required: delman, prgmods, dops */ /* end of program */ /* begin module version */ #define version 5.06 /* of rsgra.p 2007 Feb 9 2007 Feb 9, 5.06: prgmod upgrade 2007 Feb 9, 5.05: avoid rotation by 90 degrees crash because eof fails on empty wave file. 2003 Nov 18, 5.04: adjust location of figure on page. 2003 Nov 18, 5.01: label position axis 1998 Jun 19: point to rseq in see also 1993 Aug 19: minor changes 1991 May 31: rsgrap file becomes wave file and rsgrap takes range of the graph to show. origin 1986 feb 19 */ /* end module version */ /* begin module describe.rsgra */ /* name rsgra: rsequence graph synopsis rsgra(rsdata: in, picture: out, rsgrap, wave: in, marks: in, output: out) files rsdata: data file from rseq program picture: graph of rsequence in PostScript rsgrap: parameters to control the program. first line: two integers that define the from-to range to display wave: parameters to control the program. empty file means no cosine wave, otherwise the parameters of the wave are given one per line: extreme: char; h or l, the extreme high or low point on the curve defined by the wavelocation and wavebit wavelocation: real; the location in bases of the extreme wavebit: real; the location in bits of the extreme waveamplitude: real; the amplitude of the wave in bits wavelength: real; the wave length of the wave in bases marks: an empty file or a set of integers, one per line that are the locations of bases that should be specially marked on the graph. If the first line of the file begins with the letter 'b' and is followed by a real number, then this number defines the location of a bar to be placed on the graph immediately after the position given. output: messages to the user description Rsgra generates a graph of Rsequence versus position l. The user may overlay a cosine wave on the graph, to signify the orientation of the major grove, for example. See the discussion about the REMOVE feature in makelogo.p. see also rseq.p, makelogo.p author Thomas D. Schneider bugs 2007 Feb 9: Adobe Readers look at PostScript files and if the file contains a 90 rotate, the reader TURNS THE WHOLE PAGE. Stupid. To avoid this, break 90 degree rotations into two 45 degree rotations. Amazingly simple way to fool an idiotic program. But it doesn't work. So the rotation is not done now. For the same reason, the marks file cannot be empty. Use a large negative value to avoid having marks. */ /* end module describe.rsgra */ /* begin module rsgra.const */ #define nfield 4 /* size of field for printing n, the number of sites */ /* end module rsgra.const */ /* 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.59; (@ of dops.p 1990 December 5 */ /* begin module interact.const */ /* begin module string.const */ #define maxstring 2000 /* the maximum string */ /* end module string.const version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* end module interact.const version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module rsgra.filler.const */ #define fillermax 21 /* the size of the filler array for a string */ /* end module rsgra.filler.const */ /* begin module rsgra.type */ /* begin module wave.type */ /* to link several wave definitions together */ typedef struct waveparam { /* parameters to define a cosine wave */ /* define a cosine wave: */ Char extreme; /* h or l, the high or low extreme to be defined */ double wavelocation; /* the location in bases of the extreme */ double wavebit; /* the location in bits of the extreme */ double waveamplitude; /* the amplitude of the wave in bits */ double wavelength; /* the wave length of the wave in bases */ /* 2000 Nov 1. This is difficult because the length of a curve is not simply multiplied by a scale factor. BASEdashon: real; (* dashon interval in bases *) BASEdashoff: real; (* dashon interval in bases *) BASEdashoffset: real; (* dashon interval in bases *) */ double dashon; /* the size of on dashes in cm. dashon <= 0 means no dashes */ double dashoff; /* the size of on dashes in cm */ double dashoffset; /* the size of off dashes in cm */ double thickness; /* thickness of the wave in cm. <= 0 means default */ struct waveparam *next; /* the next wave */ } waveparam; /* end module wave.type version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* end module rsgra.type version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module rs.type */ typedef struct rstype { /* types of data in the data table from rseq */ long rstart, rstop; /* range of the data */ long l; /* position */ long nal, ncl, ngl, ntl; /* numbers of each base */ double rsl; /* rsequence(l) */ double rs; /* running sum of rs */ double varhnb; /* variance estimate for small sample correction */ double sumvar; /* sum of varhnb */ long nl; /* = a+c+g+t */ double ehnb; /* hg-e(n) */ Char flag; /* e = exact, a = approximate sampling statistics */ } rstype; /* end module rs.type */ /* begin module interact.type */ /* begin module string.type */ /* pointer to a string */ 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 */ Char *next; /* the next string in a series */ } string; /* end module string.type version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* end module interact.type version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module trigger.type */ typedef struct trigger { /* an object to be searched for */ string seek; /* the characters looked for */ long state; /* how close to triggering we are */ boolean skip; /* trigger not found- skip the line */ /* the trigger was found */ boolean found; } trigger; /* end module trigger.type version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module filler.type */ /* the following is an array used to fill a string. it is convenient to have it much shorter than the maxstring, so that it is easy to fill the string using procedure fillstring. the user must declare the value of constant fillermax. */ typedef Char filler[fillermax]; /* end module filler.type version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module rsgra.var */ Static _TEXT rsdata; /* output of rseq program */ Static _TEXT picture; /* output of this program, graph in PostScript */ Static _TEXT rsgrap; /* parameters to control the program */ Static _TEXT wave; /* parameters to control the drawing of the cosine wave */ Static _TEXT marks; /* locations of marks */ /* end module rsgra.var */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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 = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module interact.clearstring */ /* begin module clearstring */ /* These modules clear strings in various ways */ /* ---- */ Static Void emptystring(ribbon) string *ribbon; { /* clearstring */ /* empty the contents of the string but do NOT remove the pointer. This is useful for clearing one string within a linked list of them. */ long index; /* to the ribbon */ for (index = 0; index < maxstring; index++) ribbon->letters[index] = ' '; ribbon->length = 0; ribbon->current = 0; } /* emptystring */ /* ---- */ Static Void clearstring(ribbon) string *ribbon; { /* empty the string and remove the pointer */ emptystring(ribbon); ribbon->next = NULL; } /* clearstring */ /* ---- */ Static Void initializestring(ribbon) string *ribbon; { /* start the string with a nil pointer. This routine should be called before doing linked list work. This allows the standard string routines to clear the string without killing the pointer. This is now deprecated, do not use it since clearstring still clears the next pointer. */ printf("remove initializestring routine!\n"); printf("replace it with clearstring routine!\n"); halt(); /* to force deprecation */ clearstring(ribbon); ribbon->next = NULL; } /* initializestring */ /* end module clearstring version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* end module interact.clearstring version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module interact.writestring */ /* begin module writestring */ Static Void writestring(tofile, s) _TEXT *tofile; string *s; { /* write the string s to file tofile, no writeln */ long i; /* index to s */ long FORLIM; FORLIM = s->length; for (i = 0; i < FORLIM; i++) putc(s->letters[i], tofile->f); } /* writestring */ /* end module writestring version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* end module interact.writestring version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module trigger.proc */ /* this module allows one to scan a series of characters, as from an array or a file, and to "trigger" or detect a simple string in the series. the advantage of the trigger is that several triggers can "observe" a stream of characters at once, each looking for a different thing. some other modules required: interact.const, interact.type */ Static Void resettrigger(t) trigger *t; { /* reset the trigger to ground state */ t->state = 0; t->skip = false; t->found = false; } /* resettrigger */ Static Void testfortrigger(ch, t) Char ch; trigger *t; { /* look at the character ch. if it is part of the trigger (at the current trigger state), then the trigger state goes higher. if it is not part of the trigger then the trigger state is reset, skip is true and one should skip onward to find the trigger. if the trigger is found, found is true. 1996 Sep 12: Bug found! In the case of a trigger "ab", the program used to miss it for situations like "aab". This was because at the first a it would step up. Then it would see the second a and recognize that was not part of ab. It would fail to realize that it could be the start of a new one. The code now accounts for that possibility. */ t->state++; /* writestring(list,seek); writeln(list,'testfortrigger seek.letters[',state:1,']:', seek.letters[state],' ch:',ch); */ if (t->seek.letters[t->state - 1] == ch) { t->skip = false; if (t->state == t->seek.length) t->found = true; else t->found = false; return; } /* it failed. But wait! It could be the beginning of a NEW trigger string! */ if (t->seek.letters[0] == ch) { t->state = 1; t->skip = false; t->found = false; return; } t->state = 0; t->skip = true; t->found = false; /* reset trigger */ } /* testfortrigger */ /* end module trigger.proc version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module filler.fillstring */ Static Void fillstring(s, a) string *s; Char *a; { /* this procedure makes it reasonably easy to fill the string s with characters. one calls the procedure as: */ /* 1 2 3 4 5 */ /* 12345678901234567890123456789012345678901234567890 */ /* fillstring(s, 'this-is-the-string '); the two comments make it easy to line the characters up. also, for this example, it was assumed that the length of filler as defined by the constant fillermax was 50. */ long length = fillermax; /* of the string without trailing blanks */ long index; /* of s */ clearstring(s); while (length > 1 && a[length-1] == ' ') length--; if (length == 1 && a[length-1] == ' ') { printf("fillstring: the string is empty\n"); halt(); } for (index = 0; index < length; index++) s->letters[index] = a[index]; s->length = length; s->current = 1; } /* fillstring */ /* end module filler.fillstring version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module filler.filltrigger */ Static Void filltrigger(t, a) trigger *t; Char *a; { /* fill the trigger t */ fillstring(&t->seek, a); } /* fillstring */ /* end module filler.filltrigger version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module copyaline */ Static Void copyaline(fin, fout) _TEXT *fin, *fout; { /* copy a line from file fin to file fout */ while (!P_eoln(fin->f)) { putc(P_peek(fin->f), fout->f); getc(fin->f); } fscanf(fin->f, "%*[^\n]"); getc(fin->f); putc('\n', fout->f); } /* copyaline */ /* end module copyaline version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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 */ /* Only use initgraphics if really necessary. Since it restarts the graphics, one can't put the figure inside another PostScript figure. writeln(afile,'initgraphics'); (* 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 */ /* clean residue from before */ fprintf(afile->f, "clear erasepage %% REMOVE FOR PACKAGING INTO ANOTHER FIGURE\n"); 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.59; (@ of dops.p 1990 December 5 */ /* 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 %% REMOVE FOR PACKAGING INTO ANOTHER FIGURE\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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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: rsgra.p, line 847: * 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* ********************************************************************** */ /* end module pic.functions version = 2.59; (@ of dops.p 1990 December 5 */ /* begin module pic.degtorad */ Static double degtorad(angle) double angle; { /* convert angle in degrees to radians */ return (angle / 360 * 2 * pi); } /* end module pic.degtorad version = 2.59; (@ of dops.p 1990 December 5 */ /* 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.59; (@ of dops.p 1990 December 5 */ /* begin module pic.arc */ Static Void arc(thefile, angle1, angle2, radius, steps) _TEXT *thefile; double angle1, angle2, radius; long steps; { /* create an arc in thefile going from angle1 to angle2 (degrees) in the positive direction of angle, with the given radius. use the given number of steps to make it. return to the same position as before the arc was drawn. */ double dtheta; /* change in theta */ /* s: integer; (@ index to the steps */ double theta; /* current angle */ double x, y; /* coordinates around starting point */ /* zerox,zeroy: real; (@ starting location, center of curve */ /* zerox := picxglobal; zeroy := picyglobal; */ theta = degtorad(angle1); dtheta = degtorad(fabs(angle2 - angle1) / steps); polrec(radius, theta, &x, &y); /* can't do this for postscript arc: movea(thefile,zerox+x,zeroy+y); ' ',scale*zerox: picfield:picwidth, ' ',scale*zeroy: picfield:picwidth, */ /* NONSTANDARD postscript, much faster */ /* 'stroke newpath', */ /* force there to be no current point */ /* ' 0 0', */ fprintf(thefile->f, "a %*.*f %*.*f %*.*f\n", picfield, picwidth, scale * radius, picfield, picwidth, angle1, picfield, picwidth, angle2); fprintf(thefile->f, "arc"); if (angle2 < angle1) /* for negative draws */ putc('n', thefile->f); /* origin move: writeln(thefile,' stroke newpath 0 0 moveto'); */ fprintf(thefile->f, " n\n"); /* the moveto puts us back to the origin */ /* for s := 1 to steps do begin theta := theta + dtheta; polrec(radius,theta, x,y); linea(thefile,zerox+x,zeroy+y); end; movea(thefile,zerox,zeroy) */ } /* end module pic.arc version = 2.59; (@ of dops.p 1990 December 5 */ /* begin module pic.circler */ Static Void circler(afile, radius) _TEXT *afile; double radius; { /* make a circle at the current position of some radius. */ long steps; /* number of steps to make the circle */ /* number of segments increases with diameter, but the constant still should be a function of how good it looks on a particular graphic system, I'm afraid. However, there should be a lower bound on the number of steps, so even small circles look good */ if (radius < 1.0) steps = 25; else steps = (long)floor(radius * 25 + 0.5); arc(afile, 0.0, 360.0, radius, steps); } #define tab 9 /* tab character */ /* end module pic.circler version = 2.59; (@ of dops.p 1990 December 5 */ /* ********************************************************************** */ /* begin module skipblanks */ /* 2003 July 31: tab is considered a blank character */ Static boolean isblank(c) Char c; { /* is the character c blank or tab? */ return (c == ' ' || c == tab); } #undef tab Static Void skipblanks(thefile) _TEXT *thefile; { /* skip over blanks until a non-blank, or end of line, is found */ while (isblank(P_peek(thefile->f)) & (!P_eoln(thefile->f))) getc(thefile->f); } Static Void skipnonblanks(thefile) _TEXT *thefile; { /* skip over nonblanks until a blank, or end of line, is found */ while ((!isblank(P_peek(thefile->f))) & (!P_eoln(thefile->f))) getc(thefile->f); } Static Void skipcolumn(thefile) _TEXT *thefile; { /* skip over a data column */ skipblanks(thefile); skipnonblanks(thefile); } /* end module skipblanks version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* ********************************************************************** */ /* 2007 Feb 9 Block original routine which fails because the eof test fails (an empty afile is not detected as empty, then the read fails because it is empty). (* begin module rsgra.readwaveparameters *) procedure readwaveparameters(var afile: text; var p: waveparam); (* read from afile the parameters p, as defined by type param *) begin writeln(output,'BUBBA 0'); reset(afile); with p do begin if eof(afile) then wave := false else begin wave := true; writeln(output,'BUBBA 1'); readln(afile,extreme); readln(afile,wavelocation); readln(afile,wavebit); readln(afile,waveamplitude); readln(afile,wavelength); end end end; (* end module rsgra.readwaveparameters *) */ /* begin module readwaveparameters */ Static Void readawaveparameter(afile, wp) _TEXT *afile; waveparam *wp; { /* cmperbase, cmperbit: real); */ /* Read one wave parameter from a file. */ /* 2000 Nov 1: The following experimental code is not implemented: See the coscurve program (http://www.ccrnp.ncifcrf.gov/~toms/delila/coscurve.html) for details about the conversion from bases to cm. The code now requires cm per base and cmperbit. const lover2pi = 1.2160067233; (* L/2pi for L being the length of the cosine function from 0 to 2pi *) var multiplier: real; (* factor to multiply by to convert from bases to cm. *) */ fscanf(afile->f, "%c%*[^\n]", &wp->extreme); getc(afile->f); if (wp->extreme == '\n') wp->extreme = ' '; fscanf(afile->f, "%lg%*[^\n]", &wp->wavelocation); getc(afile->f); fscanf(afile->f, "%lg%*[^\n]", &wp->wavebit); getc(afile->f); fscanf(afile->f, "%lg%*[^\n]", &wp->waveamplitude); getc(afile->f); fscanf(afile->f, "%lg%*[^\n]", &wp->wavelength); getc(afile->f); if (wp->wavelength <= 0.0) { printf("wave parameters: wavelength must be positive\n"); halt(); } /* one may skip the dash and thickness parameters of the last wave definition. This allows backwards compatability, but should not be used in general. */ if (BUFEOF(afile->f)) { wp->dashon = 0.0; wp->dashoff = 0.0; wp->dashoffset = 0.0; wp->thickness = 0.0; return; } /* implement new method to read dash information without breaking previous wave parameter files */ if (P_peek(afile->f) != 'd') { /* default to old method */ fscanf(afile->f, "%lg%*[^\n]", &wp->dashon); getc(afile->f); wp->dashoff = wp->dashon; wp->dashoffset = 0.0; } else { getc(afile->f); /* skip the 'd' */ /* reading the parameters directly in cm would require that the user compute or guess the values. So read in as bases and compute cm. Unfortunately the length of a curve is not the simple scale factor. Make the user guess for now. */ fscanf(afile->f, "%lg%lg%lg%*[^\n]", &wp->dashon, &wp->dashoff, &wp->dashoffset); getc(afile->f); /* experimental code to account for actual wave length: readln(afile, wp.BASEdashon, wp.BASEdashoff, wp.BASEdashoffset); multiplier := 1; multiplier := lover2pi * wp.waveamplitude (* in bits *) * cmperbit (* converts to cm *) / wp.wavelength (* 10.6 bp/turn *) * cmperbase ; multiplier := lover2pi; multiplier := 1.0; multiplier := 1.0/lover2pi; multiplier := lover2pi * wp.waveamplitude (* in bits *) * cmperbit (* converts to cm *) * wp.wavelength (* 10.6 bp/turn *) * cmperbase; multiplier := 1.0 wp.dashon := multiplier * wp.BASEdashon; wp.dashoff := multiplier * wp.BASEdashoff; wp.dashoffset := multiplier * wp.BASEdashoffset; with wp do begin writeln(output, 'cmperbase = ',cmperbase:5:2); writeln(output, 'cmperbit = ',cmperbit:5:2); writeln(output, 'dashon = ',dashon:5:2); writeln(output, 'dashoff = ',dashoff:5:2); writeln(output, 'dashoffset = ',dashoffset:5:2); writeln(output, 'waveamplitude = ',wp.waveamplitude:5:2); writeln(output, 'wavelength = ',wp.wavelength:5:2); end; */ } if (BUFEOF(afile->f)) wp->thickness = 0.0; else { fscanf(afile->f, "%lg%*[^\n]", &wp->thickness); getc(afile->f); } } /* Local variables for readwaveparameters: */ struct LOC_readwaveparameters { _TEXT *afile; boolean done; /* done reading */ } ; Local Void waystoend(LINK) struct LOC_readwaveparameters *LINK; { /* if end of file, a blank line or a period, consider the job done */ boolean clear = false; /* not eof, so clear to test for comment */ if (BUFEOF(LINK->afile->f)) LINK->done = true; if (!LINK->done) { /* skip any number of comments */ while (!clear && !LINK->done) { if (BUFEOF(LINK->afile->f)) { LINK->done = true; clear = true; } else if (P_peek(LINK->afile->f) == '*') { fscanf(LINK->afile->f, "%*[^\n]"); getc(LINK->afile->f); } else clear = true; } } if (!LINK->done) { if (P_eoln(LINK->afile->f)) LINK->done = true; } if (LINK->done) return; if (P_peek(LINK->afile->f) != '.') return; fscanf(LINK->afile->f, "%*[^\n]"); getc(LINK->afile->f); LINK->done = true; } Static Void readwaveparameters(afile_, w) _TEXT *afile_; waveparam **w; { /* cmperbase, cmperbit: real); */ /* read from afile the wave parameters w. The routine is done either at end of file, when a period ('.') appears or at a blank line. This makes it compatable with different reading mechanisms. */ struct LOC_readwaveparameters V; boolean morethanone = false; /* more than one wave? */ waveparam *p; /* one of the definitions of a cosine wave */ V.afile = afile_; V.done = false; *w = NULL; waystoend(&V); if (V.done) { *w = NULL; return; } *w = (waveparam *)Malloc(sizeof(waveparam)); p = *w; while (!V.done) { waystoend(&V); if (V.done) { p->next = NULL; continue; } if (P_peek(V.afile->f) == '*') { /* skip any number of comments */ fscanf(V.afile->f, "%*[^\n]"); getc(V.afile->f); continue; } if (morethanone) { p->next = (waveparam *)Malloc(sizeof(waveparam)); p = p->next; } /* readawaveparameter(afile, p^, cmperbase, cmperbit); */ readawaveparameter(V.afile, p); /* if we ever come back again there will be more than one: */ morethanone = true; } } /* end module readwaveparameters version = 5.31; (@ of prgmod.p 2006 Oct 10 */ /* begin module rsgra.header */ Static Void header(infile, marks, outfile) _TEXT *infile, *marks, *outfile; { /* do the header of the plot */ long index; /* to lines in the infile */ if (*infile->name != '\0') { if (infile->f != NULL) infile->f = freopen(infile->name, "r", infile->f); else infile->f = fopen(infile->name, "r"); } else rewind(infile->f); if (infile->f == NULL) _EscIO2(FileNotFound, infile->name); RESETBUF(infile->f, Char); if (*marks->name != '\0') { if (marks->f != NULL) marks->f = freopen(marks->name, "r", marks->f); else marks->f = fopen(marks->name, "r"); } else rewind(marks->f); if (marks->f == NULL) _EscIO2(FileNotFound, marks->name); RESETBUF(marks->f, Char); if (*outfile->name != '\0') { if (outfile->f != NULL) outfile->f = freopen(outfile->name, "w", outfile->f); else outfile->f = fopen(outfile->name, "w"); } else { if (outfile->f != NULL) rewind(outfile->f); else outfile->f = tmpfile(); } if (outfile->f == NULL) _EscIO2(FileNotFound, outfile->name); SETUPBUF(outfile->f, Char); /* writeln(outfile,'.nf'); writeln(outfile,'% rsgra ',version:4:2); */ fprintf(outfile->f, "%%!PS-Adobe-2.0 EPSF-2.0\n"); fprintf(outfile->f, "%%%%Title: rsgra %4.2f\n", version); fprintf(outfile->f, "%%%%Creator: Tom Schneider\n"); fprintf(outfile->f, "%%%%BoundingBox: 0 0 630 900\n"); fprintf(outfile->f, "%%%%DocumentFonts:\n"); fprintf(outfile->f, "%%%%EndComments\n"); fprintf(outfile->f, "%%%%EndProlog\n"); /* copy header lines to outfile */ for (index = 1; index <= 3; index++) { fprintf(outfile->f, "%% "); copyaline(infile, outfile); } fprintf(outfile->f, "\n%% MARKS: \n"); /* identify the marks */ while (!BUFEOF(marks->f)) { fprintf(outfile->f, "%% "); copyaline(marks, outfile); } if (*marks->name != '\0') { if (marks->f != NULL) marks->f = freopen(marks->name, "r", marks->f); else marks->f = fopen(marks->name, "r"); } else rewind(marks->f); if (marks->f == NULL) _EscIO2(FileNotFound, marks->name); RESETBUF(marks->f, Char); } #define step 0.01 /* step size for the lines drawn */ /* Local variables for cosine: */ struct LOC_cosine { double amplitude, phase, wavelength, base; } ; Local double c(y, LINK) double y; struct LOC_cosine *LINK; { /* do the actual calculation of the cosine given y */ return (LINK->amplitude * cos((-y - LINK->phase) / LINK->wavelength * 2 * pi) + LINK->base); } Local boolean between(a, b, c, LINK) double a, b, c; struct LOC_cosine *LINK; { /* is b between a and c? */ return (a <= b && b <= c); } /* end module rsgra.header */ /* begin module rsgra.cosine */ Static Void cosine(afile, amplitude_, phase_, wavelength_, base_, xmin, ymin, xmax, ymax) _TEXT *afile; double amplitude_, phase_, wavelength_, base_, xmin, ymin, xmax, ymax; { /* draw a cosine function down the page: x := amplitude*cos( ((-y-phase)/wavelength)*2*pi) + base, where the given parameters are in inches. parts outside the box defined by xmin,ymin,xmax,ymax are not drawn. pi is defined in pic.const */ struct LOC_cosine V; boolean inside1; /* point 1 was inside the bounds */ boolean inside2; /* point 2 was inside the bounds */ double x1, y1, x2, y2; /* internal coordinates */ V.amplitude = amplitude_; V.phase = phase_; V.wavelength = wavelength_; V.base = base_; /* start graphing at the top of the plot */ y1 = ymax; x1 = c(y1, &V); inside1 = between(ymin, y1, ymax, &V) & between(xmin, x1, xmax, &V); while (y1 >= ymin) { x2 = x1; y2 = y1; inside2 = inside1; y1 -= step; x1 = c(y1, &V); inside1 = between(ymin, y1, ymax, &V) & between(xmin, x1, xmax, &V); if (inside1 && inside2) { movea(afile, x1, y1); linea(afile, x2, y2); } } } #undef step /* end module rsgra.cosine */ /* ********************************************************************** */ /* ********************************************************************** */ /* begin module rs.readrsrange */ Static Void readrsrange(rsdata, r) _TEXT *rsdata; rstype *r; { /* read the range data from rsdata to r. data is assumed to be the rsdata file from program rseq. */ long index; /* for counting lines of rsdata */ Char skip; /* a character to skip the '*' on the line */ for (index = 1; index <= 11; index++) { fscanf(rsdata->f, "%*[^\n]"); getc(rsdata->f); } fscanf(rsdata->f, "%c%ld%ld%*[^\n]", &skip, &r->rstart, &r->rstop); getc(rsdata->f); /* writeln(output, 'range: ',r.rstart:1,' ',r.rstop:1);*/ if (skip == '\n') skip = ' '; } /* end module rs.readrsrange */ /* begin module rs.getrsbegin */ Static Void getrsbegin(infile) _TEXT *infile; { /* skip to the beginning of the data in a data file from the rseq program */ Char ch; /* a character read from infile */ trigger begindata; /* a trigger to locate the beginning of the data */ /* 1 2 */ /* 123456789012345678901 */ filltrigger(&begindata, "l a c g t"); resettrigger(&begindata); if (*infile->name != '\0') { if (infile->f != NULL) infile->f = freopen(infile->name, "r", infile->f); else infile->f = fopen(infile->name, "r"); } else rewind(infile->f); if (infile->f == NULL) _EscIO2(FileNotFound, infile->name); RESETBUF(infile->f, Char); while (!begindata.found) { if (P_eoln(infile->f)) { fscanf(infile->f, "%*[^\n]"); getc(infile->f); } if (BUFEOF(infile->f)) { printf("beginning of data not found\n"); halt(); } ch = getc(infile->f); if (ch == '\n') ch = ' '; testfortrigger(ch, &begindata); } fscanf(infile->f, "%*[^\n]"); getc(infile->f); } /* end module rs.getrsbegin */ /* begin module rs.readrsdata */ Static Void readrsdata(rsdata, rdata) _TEXT *rsdata; rstype *rdata; { /* read data from the data file of program rseq into the datatype */ fscanf(rsdata->f, "%ld%ld%ld%ld%ld%lg%lg%lg%lg%ld%lg", &rdata->l, &rdata->nal, &rdata->ncl, &rdata->ngl, &rdata->ntl, &rdata->rsl, &rdata->rs, &rdata->varhnb, &rdata->sumvar, &rdata->nl, &rdata->ehnb); /* skip spaces to find the flag: */ while (P_peek(rsdata->f) == ' ') getc(rsdata->f); fscanf(rsdata->f, "%c%*[^\n]", &rdata->flag); getc(rsdata->f); /* writeln(output,'readrsdata: l a c g t flag = ', l:1,' ',nal:1,' ',ncl:1,' ',ngl:1,' ',ntl:1,' ',flag); */ if (rdata->flag == '\n') rdata->flag = ' '; } #define duplicates 2 /* number of times to draw the main line */ #define halfdatawidth 2.4 /* inches to shift the data to the left from the left edge of the graph. the lowpoint is added to this to find the dataposition */ #define leftbound 3.0 /* distance of the left bound of the graph from the left side of the paper, (not including the data) */ #define spacing 0.01 /* spacing between the multiple drawings of the main line */ #define xscale 2.00 /* inches per bit across the plot */ #define yscale 0.1595 /* end module rs.readrsdata */ /* ********************************************************************** */ /* ********************************************************************** */ /* begin module themain */ Static Void themain(infile, rsgrap, wave, marks, outfile) _TEXT *infile, *rsgrap, *wave, *marks, *outfile; { /* the main procedure of the program */ /* vertical spacing between graph lines, ie, inches per base */ /* 0.15 inches worked nicely for a while, until I wanted the graph to match the standard printer output (unix pr). There were 39 bases in the range, and they should have been 2.3 bases longer, so the new yscale is 0.15 * (39+2.3)/39. */ double bar; /* the size of an error bar */ long barposition; /* position after which to write a bar */ long basetomark; /* the base to be marked according to the marks file */ double dataposition; /* amount to shift the data to the left from the zero of the graph */ double dx, dy; /* changes of x and y to make multiple lines */ long frombase, tobase; /* range of the graph to display */ long line; /* an index to the multiple lines drawn */ double linedx, linedy; /* line * dx or dy, the actual jump made for multiple lines (done for speed) */ double lowpoint = LONG_MAX; /* the lowest point in the entire plot */ boolean notfirstpoint = false; /* true for all points except the first point plotted */ /* params: waveparam; (* parameters to control the program *) */ waveparam *params; /* wave parameters to control the program */ double point; /* a low point in the graph */ long position; /* a location in the aligned sequence */ rstype rdata; /* data from rseq */ double xmin, ymin, xmax = 2.0 * xscale, ymax = 0.0; /* the bounds of the graph, in inches */ double x1, y1; /* the current end of the curve */ double x2, y2; /* the previous end of the curve */ double ybar; /* location of the bar in graphics coordinates */ long FORLIM; printf("rsgra %4.2f\n", version); /* copy header stuff */ header(infile, marks, outfile); if (*wave->name != '\0') { if (wave->f != NULL) wave->f = freopen(wave->name, "r", wave->f); else wave->f = fopen(wave->name, "r"); } else rewind(wave->f); if (wave->f == NULL) _EscIO2(FileNotFound, wave->name); RESETBUF(wave->f, Char); readwaveparameters(wave, ¶ms); if (*rsgrap->name != '\0') { if (rsgrap->f != NULL) rsgrap->f = freopen(rsgrap->name, "r", rsgrap->f); else rsgrap->f = fopen(rsgrap->name, "r"); } else rewind(rsgrap->f); if (rsgrap->f == NULL) _EscIO2(FileNotFound, rsgrap->name); RESETBUF(rsgrap->f, Char); if (BUFEOF(rsgrap->f)) { frombase = -LONG_MAX; tobase = LONG_MAX; } else { fscanf(rsgrap->f, "%ld%ld%*[^\n]", &frombase, &tobase); getc(rsgrap->f); } /* now read the rest of the parameters */ if (BUFEOF(marks->f)) barposition = -LONG_MAX; /* don't do the bar */ else { if (P_peek(marks->f) == 'b') { getc(marks->f); /* skip the 'b' */ fscanf(marks->f, "%ld%*[^\n]", &barposition); getc(marks->f); printf("bar placed at position %ld\n", barposition); } } /* get the first base to mark */ if (!BUFEOF(marks->f)) { fscanf(marks->f, "%ld%*[^\n]", &basetomark); getc(marks->f); } else basetomark = -LONG_MAX; if (*infile->name != '\0') { if (infile->f != NULL) infile->f = freopen(infile->name, "r", infile->f); else infile->f = fopen(infile->name, "r"); } else { /* find the range of the graph in bases */ rewind(infile->f); } if (infile->f == NULL) _EscIO2(FileNotFound, infile->name); RESETBUF(infile->f, Char); readrsrange(infile, &rdata); /* if params.wave */ if (params != NULL) fprintf(outfile->f, "%% cosine wavelength is %4.2f bases\n", params->wavelength); if (frombase < rdata.rstart) frombase = rdata.rstart; if (tobase > rdata.rstop) tobase = rdata.rstop; /* find lowest bound of the graph, in bits */ getrsbegin(infile); FORLIM = rdata.rstop; for (position = rdata.rstart; position <= FORLIM; position++) { /* skip lines with an '*' */ if (P_peek(infile->f) != '*') { readrsdata(infile, &rdata); point = rdata.rsl - sqrt(rdata.varhnb); if (point < lowpoint) lowpoint = point; } else { fscanf(infile->f, "%*[^\n]"); getc(infile->f); } } /* make lowpoint be just below zero, or the next tenth of a bit lower: */ if (lowpoint >= 0.0) lowpoint = -0.1; else lowpoint = (long)(1 - 10 * lowpoint) / -10.0; dataposition = halfdatawidth - lowpoint * xscale; /* define the bounds of the graph region, in inches */ xmin = lowpoint * xscale; ymin = (frombase - tobase - 2) * yscale; /* startpic(outfile,defscale,3.3,8.5,'c'); */ startpic(outfile, (double)defscale, 1.0, 5.0, 'c'); fprintf(outfile->f, "45 rotate\n"); fprintf(outfile->f, "45 rotate\n"); /* the coordinate (0,0) is the upper left hand corner of the graph, not including the data. First create the location of the data: */ movea(outfile, -leftbound, 0.0); movea(outfile, -dataposition, 0.0); fprintf(outfile->f, "gsave 1 2 scale\n"); fprintf(outfile->f, "( l a c g t Rs(l)) x\n"); fprintf(outfile->f, "grestore\n"); /* placement of the label in bits along the Y axis: */ movea(outfile, 0.25 * xscale, 3 * yscale); fprintf(outfile->f, "gsave 1 2 scale\n"); fprintf(outfile->f, "(Rs(l) = Rsequence(l), Information in bits) x\n"); fprintf(outfile->f, "grestore\n"); /* the center */ /* placement of the label in bits along the X axis: */ /* movea(outfile,-1.49*xscale,-8.0*yscale); zzz */ /* adjust for size of the words, in cm */ movea(outfile, -1.60 * xscale, ((frombase - tobase) / 2.0 + 8.5) * yscale); fprintf(outfile->f, "gsave 2 2 scale\n"); /* 2007 Feb 9: Adobe Readers look at PostScript files and if the file contains a 90 rotate, the reader TURNS THE WHOLE PAGE. Stupid. To avoid this, break 90 degree rotations into two 45 degree rotations. Amazingly simple way to fool an idiotic program. writeln(outfile,'-90 rotate'); */ /* that didn't work. Don't rotate figure writeln(outfile,'-45 rotate'); writeln(outfile,'-45 rotate'); */ fprintf(outfile->f, "-45 rotate\n"); fprintf(outfile->f, "-45 rotate\n"); fprintf(outfile->f, "(Position L (in bases)) x\n"); fprintf(outfile->f, "grestore\n"); /* now put the book id onto the graph */ /* x coordinate was xscale-dataposition, but this way more title will fit on: */ movea(outfile, -2 * dataposition / 2, 4 * yscale); putc('(', outfile->f); fprintf(outfile->f, "rsgra %4.2f ", version); if (*infile->name != '\0') { if (infile->f != NULL) infile->f = freopen(infile->name, "r", infile->f); else infile->f = fopen(infile->name, "r"); } else rewind(infile->f); if (infile->f == NULL) _EscIO2(FileNotFound, infile->name); RESETBUF(infile->f, Char); fscanf(infile->f, "%*[^\n]"); getc(infile->f); fscanf(infile->f, "%*[^\n]"); getc(infile->f); while (!P_eoln(infile->f)) { /* like copyaline but no writeln */ putc(P_peek(infile->f), outfile->f); getc(infile->f); } fprintf(outfile->f, ") x %% NOHEADER FOR PACKAGING INTO ANOTHER FIGURE\n"); /* create upper x axis */ movea(outfile, 0.0, 0.0); xaxis(outfile, xmax, 0.0, 1.0, 2.0, -0.1, 0.0, 0.1, 3L, 1L); /* more ticks: */ mover(outfile, lowpoint * xscale, 0.0); xaxis(outfile, (2.0 - lowpoint) * xscale, lowpoint, 0.1, 2.0, -0.1, 0.0, 0.0, 0L, 1L); /* dashed line down the plot: */ movea(outfile, 0.0, 0.0); drawa(outfile, 0.0, ymin, '-', 0.03); /* create lower x axis */ xaxis(outfile, xmax, 0.0, 1.0, 2.0, 0.1, 0.0, -0.1, 0L, 1L); /* more ticks: */ mover(outfile, xmin, 0.0); xaxis(outfile, (2.0 - lowpoint) * xscale, lowpoint, 0.1, 2.0, 0.1, 0.0, 0.0, 0L, 1L); /* xaxis(outfile, axlength,fromtic,interval,totic: real; length, dx, dy: real; width, decimal: integer); */ /* about to draw the first point, so there is no previous point */ /* prepare for reading the data for graphing */ getrsbegin(infile); FORLIM = rdata.rstop; /* create the curve of rsequence(l) */ for (position = rdata.rstart; position <= FORLIM; position++) { /* skip lines with an '*' */ if (P_peek(infile->f) != '*') { readrsdata(infile, &rdata); if (position >= frombase && position <= tobase) { /* move to the next data line and draw a bar */ if (position == barposition) { /* the factor multiplying dataposition determines the left end of the bar. It is made to fit nicely. */ /* the extra 0.5 makes sure that the position is between bases, so it won't overwrite the marks */ ybar = (frombase - position - 1.5) * yscale; /* movea(outfile,-1.60*dataposition,ybar); */ movea(outfile, -dataposition, ybar); linea(outfile, 2 * xscale, ybar); } /* move to the next data line and give the data */ /* the extra amount added to the position places the numbers in a better position */ movea(outfile, -dataposition, (frombase - position - 1.2) * yscale); fprintf(outfile->f, "(%*ld %*ld %*ld %*ld %*ld%7.2f) x\n", nfield, rdata.l, nfield, rdata.nal, nfield, rdata.ncl, nfield, rdata.ngl, nfield, rdata.ntl, rdata.rsl); if (notfirstpoint) { /* multiple draw line from previous point */ x2 = x1; /* set (x2,y2) to previous point */ y2 = y1; x1 = rdata.rsl * xscale; /* update (x1,y1) */ y1 = (frombase - position - 1) * yscale; /* the slope of the line from (x2,y2) to (x1,y1) is m := (y1-y2)/(x1-x2); the multiple lines must be drawn perpendicular to this, m' := -1/m; the arctan of this gives the angle to move the lines in. conversion from polar to rectangular coordinates completes the job. note that y2-y1 will never be zero. */ polrec(spacing, atan((x1 - x2) / (y2 - y1)), &dx, &dy); for (line = -duplicates; line <= duplicates; line++) { linedx = line * dx; linedy = line * dy; movea(outfile, x2 + linedx, y2 + linedy); linea(outfile, x1 + linedx, y1 + linedy); } } else { x1 = rdata.rsl * xscale; y1 = (frombase - position - 1) * yscale; notfirstpoint = true; } /* the first point does not have a line drawn to it */ /* create error bar and circle */ /* move back to the location of the new data point: */ movea(outfile, x1, y1); bar = sqrt(rdata.varhnb) * xscale; mover(outfile, bar, 0.0); liner(outfile, -bar, 0.0); circler(outfile, yscale / 3); liner(outfile, -bar, 0.0); /* create specially marked bases */ if (basetomark == rdata.l) { /*movea(outfile,lowpoint+xscale*0.025,*/ movea(outfile, lowpoint, (frombase - position - 1) * yscale); circler(outfile, yscale / 3); printf("marking base %ld\n", rdata.l); if (!BUFEOF(marks->f)) { fscanf(marks->f, "%ld%*[^\n]", &basetomark); getc(marks->f); } else basetomark = -LONG_MAX; } } } else { fscanf(infile->f, "%*[^\n]"); getc(infile->f); } } /* overlay a cosine wave as specified by the parameters */ /* if params.wave then with params do begin */ if (params != NULL) { if (params->extreme == 'h') /* high center point cosine wave */ cosine(outfile, xscale * params->waveamplitude, yscale * (params->wavelocation - frombase + 1), yscale * params->wavelength, xscale * (params->wavebit - params->waveamplitude), xmin, ymin, xmax, ymax); else cosine(outfile, -xscale * params->waveamplitude, yscale * (params->wavelocation - frombase + 1), yscale * params->wavelength, xscale * (params->wavebit + params->waveamplitude), xmin, ymin, xmax, ymax); } /* low center point cosine wave */ stoppic(outfile); } #undef duplicates #undef halfdatawidth #undef leftbound #undef spacing #undef xscale #undef yscale /* end module themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; marks.f = NULL; strcpy(marks.name, "marks"); wave.f = NULL; strcpy(wave.name, "wave"); rsgrap.f = NULL; strcpy(rsgrap.name, "rsgrap"); picture.f = NULL; strcpy(picture.name, "picture"); rsdata.f = NULL; strcpy(rsdata.name, "rsdata"); themain(&rsdata, &rsgrap, &wave, &marks, &picture); _L1: if (rsdata.f != NULL) fclose(rsdata.f); if (picture.f != NULL) fclose(picture.f); if (rsgrap.f != NULL) fclose(rsgrap.f); if (wave.f != NULL) fclose(wave.f); if (marks.f != NULL) fclose(marks.f); exit(EXIT_SUCCESS); } /* End. */