/code_comm/RecentFuncs.c

/* tacg - a command line tool for the analysis of nucleic acids and protein  */

/* Copyright Š 1994-2001 Harry J Mangalam, tacg Informatics

(mangalam@tacgi.com | mangalam@home.com,  | 949 856 2847) */



/* $Id: RecentFuncs.c,v 1.18 2001/11/05 02:48:42 hjm Exp $  */



/* The use of this software (except that by Harald T. Alvestrand, which is described in 'udping.c')

   and that by James Knight, which is described in 'seqio.c' is bound by the notice that appears 

   in the file 'tacg.h' which should accompany this file.  In the event that 'tacg.h' is not bundled 

   with this file, please contact the author.

*/



#include <stdio.h>

#include <ctype.h> 

#include <string.h>

#include <stdlib.h>

#include <regex.h>

#include <math.h>

#include <sys/types.h>

#include <sys/stat.h>

#include <unistd.h>



#include "seqio.h"

#include "tacg.h"  

/* #include <pcre.h>  */



/* GrafMap is the 1st take on the long-delayed graphical plasmid/linear maps function */

/* on calling, runs thru RE and makes a 1st pass guestimate of how to map all the hits. If there

	are too many hits, it should truncate at a level which does not leave a pattern partially 

   mapped.   As can be determined by examining the Postscript generated, I'm a PS newbie and the PS 

   could be made significantly more compact by encoding the tic routines in PS for example and letting

   the PS engine do all the calculations. It's also full of debug info, which will remain until the

   code solidifies more.  The PDF output tho is quite robust.

*/

   

void GrafMap(long *Degen_Log, int NumProtos, int *Protos, char *Description, char datestamp[80], float FrACGT[4]) {



   int i, j, n, bi, ti, sli, HPc, REi, slotted, FCnt, ORF_line_width, f, ORFNum, 

   	skip_this_map=0, need_prolog=0, NLabels, RWCincr, previ, start, end, run, 

      sep, center, abegin, HPEnd = NUM_PLASMID_LABELS+30,

       

      debug=0, dstart=0, dend=15; /* set the debug vars */

       

   long HPi=0, HPMax=0, RWC, ldiff, BasesPerLabel,

        seq_len = F.Seq_Len, temp, 

        HP[3][HPEnd], /* Holding Pen for counting sites 0=site 1 = REi 2=adj*/

        mdbl; /*  min dist between labels */ 

   char  *MapFileName,

         *ptmp,

         *rstr, *o, star, *arrow, 

         Label[20];

   float ORF_density, x_left, x_right, y_top, y_bot, yoffset;

   /* some constants, measurements in inches, arcs in degrees, */

   double Radius 		= 3.0, 

		   CenterX 		= 4.25, 

         CenterY 		= 4.5,

         /* xold, yold, */

         AnglePerLabel, 

         /* pi = 3.141592653589, */

         topx, topy, arc_end, hairline_radius,

         PtCenterX, PtCenterY, PtRadius;

         

   FILE *PS;     /* output postscript file */

   FILE *PROLOG; /* prolog file */



   struct NextHit_struct NH; 

   struct stat *Fbuf;

   

   /* mangle the file name to be compatible with a cgi call if it was so called */

   if ((MapFileName = calloc(64, sizeof(char))) == NULL) { BadMem("RecentFuncs:MapFileName", 1); }

   if (F.HTML != -1) { /* then we need the /tmp-named files */

      sprintf(MapFileName, "%s/tacg_Map.ps", F.tmppath);

   } else {

      MapFileName = "tacg_Map.ps"; /* name by itself to write into the $PWD */

   }

   PtCenterX = CenterX * 72; 

   PtCenterY = CenterY * 72; 

   PtRadius  = Radius  * 72;

   x_left = y_bot = 25; /* ~1/3" up from the margin to start */

   x_right = 576; /* ~8" right margin */

   y_top = 756; /* ~10.5" top margin */

   

   NLabels = NUM_PLASMID_LABELS; /*  how many Label slots there are */

   mdbl = seq_len / (long)NLabels;

   AnglePerLabel = 360 / NUM_PLASMID_LABELS;

   BasesPerLabel = (long)((long)AnglePerLabel * seq_len / 360);

  

   for (i=0;i<HPEnd;i++) { for (j=0; j<3; j++) {HP[j][i] = 0;}  }

   

	if ((rstr = calloc(MAX_LEN_RULESTRING,sizeof(char))) == NULL) { BadMem("Can't init mem for 'rstr'", 1); }

   if ((Fbuf = calloc(1, sizeof(*Fbuf))) == NULL) { BadMem("RecentFuncs:Fbuf",1); }

   

   if ((ptmp = (SearchPaths("psprolog", "PSPROLOG file"))) == NULL) {

      fprintf(stderr,"Can't open the 'psprolog' file!! - Check spelling, ownership, existance, etc - Bye!!\n");

      exit (1);

   } 

   if ((PROLOG=fopen(ptmp,"r")) == NULL) {  /* or the standard/optional one from SetFlags() */

      fprintf(stderr,"Cannot open the postscript prolog file 'psprolog' for reading!!\n"); 

      exit(1); /* print an error and die gracefully */

   }      

      

   /* if file doesn't exist, need to write in prolog */

   if (stat(MapFileName, Fbuf) == -1) { need_prolog = 1; }

   

   if ((PS=fopen(MapFileName,"a")) == NULL) {  /* or the standard/optional one from SetFlags() */

      fprintf(stderr,"Cannot open the output postscript file \"%s\" for writing!!\n", MapFileName); 

      exit(1); /* print an error and die gracefully */

   } 

   

   if (need_prolog) {

      /* use the Illustrator PS prolog as it has a number of useful defs, so copy PROLOG into PS */

      while (fgets(rstr, 512, PROLOG)) { fprintf(PS, "%s", rstr);  }

   }

   

   /* ------------------------  write the title, date, etc --------------------------------------- */

   /* write a header at top left */

   yoffset = y_top;

   fprintf(PS, "\n\n/Helvetica-Bold findfont 12 scalefont setfont\n");

   fprintf(PS, "newpath %f %f moveto\n (Description: %s) show      \n", 

                     x_left, yoffset, Description); 

   yoffset -= 12;

   fprintf(PS, "%f %f moveto\n (Date: %s) show\n", x_left, yoffset, datestamp );  

                

   fprintf(PS, "\n\n/Helvetica findfont 10 scalefont setfont\n");

   yoffset -= 10;

   fprintf(PS, "newpath %f %f moveto\n (Map by tacg v. %s (http://24.0.123.214/tacg3)) show\n", 

                     x_left, yoffset, TACG_VERSION);

   fprintf(PS, "\n\n/Helvetica findfont 10 scalefont setfont\n");

   yoffset -= 10;

   fprintf(PS, "newpath %f %f moveto\n (Suggestions to Harry Mangalam (mangalam@home.com)) show\n", 

                     x_left, yoffset);

    

   /* print a little legend at bottom left */

   yoffset = 30;

   fprintf(PS, "\n/Helvetica findfont 8 scalefont setfont\n");

   fprintf(PS, "newpath %f %f moveto\n (Label = Name[base#]overhang (*), where: ) show\n", x_left, (y_bot+yoffset));

   yoffset -= 8;

   fprintf(PS, " %f %f moveto\n (Name = Name of the RE or pattern,  ) show\n", x_left, (y_bot+yoffset));

   yoffset -= 8;

   fprintf(PS, " %f %f moveto\n ([base#] = cut site of the RE or ~center of the pattern) show\n",  x_left, (y_bot+yoffset));

   yoffset -= 8;

   fprintf(PS, " %f %f moveto\n (overhang: \\\\ = 5', / = 3', | = blunt. '*' = 1 hit only.) show\n",  x_left, (y_bot+yoffset));



   /* assume an 8x10.5 US letter page; will probably have to correct to various media later */

   

   /* set the line width and then draw a circle at the center of the page */

   fprintf(PS, "1 setlinewidth\nnewpath %f %f %f 0 360 arc stroke\n", PtCenterX, PtCenterY, PtRadius);

   /* outline circle to set off degens */

   fprintf(PS, "0.3 setlinewidth\nnewpath %f %f %f 0 360 arc stroke\n", PtCenterX, PtCenterY, PtRadius+3);

   

   /* write a small + at the center for a registration mark */

   fprintf(PS, " .3 setlinewidth %f %f moveto\n 10 0 rlineto stroke\n %f %f moveto\n 0 10 rlineto stroke\n", 

      (PtCenterX - 5), PtCenterY, PtCenterX, (PtCenterY - 5));

	



/* --------------------------   and the sequence stats -----------------------------------  */



   fprintf(PS, "\n/Helvetica-Bold findfont 10 scalefont setfont\n");   

   fprintf(PS, "newpath %f %f moveto\n", (PtCenterX - 30), PtCenterY+30);

   fprintf(PS, " (%ld bases)  show\n", seq_len) ;





   fprintf(PS, "\n/Helvetica findfont 6 scalefont setfont\n");   

   fprintf(PS, "newpath %f %f moveto\n", (PtCenterX - 100), PtCenterY+10);

   fprintf(PS, " (%ld A(%.2f %%)  %ld C(%.2f %%)  %ld G(%.2f %%)  %ld T(%.2f %%))  show\n", 

		      	F.NBases[0], FrACGT[0]*100, F.NBases[1], FrACGT[1]*100, F.NBases[2], 

               FrACGT[2]*100, F.NBases[3], FrACGT[3]*100);



   

if (seq_len > (temp = F.NBases[0] + F.NBases[1] +F.NBases[2] + F.NBases[3])) {

	fprintf(PS, "newpath %f %f moveto\n", (PtCenterX - 90), PtCenterY-15);

   fprintf (PS, "\n(NB: sequence length > A+C+G+T due to %ld IUPAC degeneracies.) show\n", seq_len-temp);

   fprintf(PS, "newpath %f %f moveto\n", (PtCenterX - 80), PtCenterY-22);

   fprintf (PS, "(# of:  N:%ld  Y:%ld  R:%ld  W:%ld  S:%ld  K:%ld  M:%ld  B:%ld  D:%ld  H:%ld  V:%ld) show \n",

        F.NBases[14], F.NBases[5], F.NBases[6], F.NBases[7], F.NBases[8], F.NBases[9], F.NBases[10], 

        F.NBases[10], F.NBases[11], F.NBases[12], F.NBases[13]);

}

     

   /* write out the alternating light/dark 1KB bands ?  not yet.. */



   /* start counting the number of sites and enter them into a local array */

   HPc = 0; j=2;

   RWCincr = seq_len / NLabels;

   while (HPc < NLabels && j < (*DD).Dat[0] && ((*DD).Dat[j] < seq_len)) {

      RWC = (*DD).Dat[j];

      REi = (int)labs((*DD).Dat[j+1]);

      if (RE_Filters_OK(REi)) { /* if the RE is allowed by the filter */

/*          fprintf(stderr,"%s accepted!\n", RE[(int)labs((*DD).Dat[j+1])].E_nam); */

      /* try to place the info in the appro slot in HP */

         slotted = 0;

                               

            HPi = (int)((float)(*DD).Dat[j] / (float)seq_len * NLabels); /* HPi points to the relative position of the hit */

            /* bc circular sequences are padded, possible to find hits 

            	before and beyond actual sequence, so have to check for over/underruns */

/*          if (HPi > HPEnd) { HPi = HPEnd;}

            if (HPi < 0) { HPi = 0;}  

*/

            if (HP[0][HPi] == 0) {     /* if there isn't anything already there */

               HP[0][HPi] = RWC; /* the RWC of the hit */ 

               HP[1][HPi] = REi;          /* the RE index */

               ldiff = HP[0][HPi-1] - HP[0][HPi];

               if (HPc > 0 && ldiff > 0 && ldiff < mdbl) {

                  HP[2][HPi] += RWCincr;

               }

               slotted = 1;

            } else { /* have to find a nearby open slot in the appro position */

               /* since this is taken from DD.Dat, the sites WILL BE in order, so should just 

                  have to step to the NEXT open slot */

      /* this routine does not guarantee that the lines won't cross, altho it should work in most cases 

         in order to guarantee this, would have to make it into a linked list structure - possible

         but lets wait to see how it behaves in real life. */

/*                fprintf(stderr,"\nCOLLISION @ %d!\n", HPi); */

               while (HP[0][HPi] != 0 && HPi < HPEnd) { /* incr until it finds an empty slot */

                  HPi++;

               }

               

               if (HPi >= HPEnd) {

               	/* then we ran into the end, but there are probably empty slots left before this

                  	so let's see if we can find them */

                  sli=HPi;

                  while (HP[0][HPi] != 0) {sli--;}

                  if (sli <= 0) {   /* we really did use them all up, so die */

                     fprintf(stderr, "\nLooks like we ran out of Label Slots in GrafMap:HP\n");

                     skip_this_map = 1;

                     exit(1); /* should skip, not exit, to handle multiple maps, but not yet */

               	} else { /* we found one before [0], but to use it, have to move all the 

                  	          trailing data UP 1 so move everything from here to end up one */

                     

                     ti=sli; bi=sli+1;

                     while (bi < HPEnd) {

                     	for (i=0; i<3; i++) { HP[i][ti++] = HP[i][bi++]; }

                     }

                     /* so we should have moved everything up by 1..? 

                     and HPi should still be pointing at the same (previously occupied) el. */

                  }

               }

               HP[0][HPi] = RWC;  /* the RWC of the hit */ 

               HP[1][HPi] = REi;  /* the RE index */  

               HP[2][HPi] += RWCincr; /* the increment caused by the too-close labels */

               slotted = 1;    

/*                fprintf(stderr,"\nHPi relocated to %d \n", HPi);          */

            }

         HP[2][HPi] += HP[0][HPi]; /*  this makes the [2] teh entire diff instead of the increment */

         HPc++;

      } else {

/*          fprintf(stderr,"\t%s rejected!\n", RE[(int)labs((*DD).Dat[j+1])].E_nam); */

      }

      j+=2;

   }

   /* this should end with the (size NLabels) Holding Pen (HP) populated with 

      RWC's (HP[0][x]) and REi's (HP[1][x]) of REs / patterns */

   

   if (HPi == HPEnd ) {

      fprintf(stdout, "\n\tMore than %d patterns met crits for plasmid map.\n"

       "\tIncrease the stringency and try again.\n", NLabels);

      exit(1);

   } else { HPMax = HPi + 1;}



if (debug > 0) {

/*  --------------------  DEBUG----------------------   */

   fprintf(stderr, "\n\nPass %d\n", debug++);

   for (i=dstart; i<dend; i++) {  

      fprintf(stderr, "\n%3d: ", i);

      for (j=0; j<3; j++) { fprintf(stderr, " %5ld ", HP[j][i]); }

   }

/*       -----------------------DEBUG -------------------------*/

}



   previ = 0;

   /* 2nd pass - increase spacing of labels that have mapped too closely  */

	HPi = 1;

   while (HPi < HPEnd) {

/*    for (HPi=1; HPi<HPEnd; HPi++)  */   /* original */

      if ( HP[0][HPi] != 0 && HP[0][HPi+1] != 0) { /* skip if not 2 in a row */

         start=HPi;

         while ((HPi < HPEnd)  && ((HP[0][HPi] != 0) || ((HP[0][HPi] != 0) && (HP[0][HPi+1] - HP[0][HPi] > mdbl)))) {



         /* could also examine the run itself to detect those spots where instead of taking the

         	next position, they are remapped to a spot corresponding to their original offset -

            this happens at the end of a long run where the last few els are actually separated 

            by a long distance  from the rest of the run, but get mapped into the run bc they

            aren't separated by any empty els */

            

/*             fprintf(stderr, "\nHPi = %ld is in a run \n", HPi); */

           previ++;

           HPi++; 

         } /* how many in a row? */

         end = --HPi; run = end - start + 1; 

         

			  if (debug > 0) {

           		fprintf(stderr, "\n at HPi= %ld, start=%d end=%d run = %d\n",   HPi, start, end, run ); 

           		for (j=0; j<3; j++) { fprintf(stderr, " %5ld ", HP[j][HPi]); }

           }



         /* calculate the separation betw labels to be equal */

         sep = (HP[2][end] - HP[2][start] ) / run;  /* sep = average separation in bp distance */

         center = ((HP[0][end] - HP[0][start]) / 2) + HP[0][start];

        abegin = center - (BasesPerLabel * (run/2)) ; /* back up the right amount */

        if (run == 1) {abegin = abegin - (BasesPerLabel*0.5);}

        if (run > 2) {abegin -= BasesPerLabel;} /* if have a larger run, back up some more */

        if (abegin < 0) { abegin = 0; }

        HP[2][start] = abegin;  /* + BasesPerLabel; mark the 1st at the backup position */

      /*  abegin += BasesPerLabel;  was sep */

        for (i = start+1; i <= end; i++) { 

            abegin += BasesPerLabel; 

            HP[2][i] = abegin; /* they should now be spread evenly across the group spread */

         }  /* and the NEXT round should detect a collision with this one. */ 

         previ = HPi;  

      }

      HPi++;

   } 



if (debug > 0) {

/*  --------------------  DEBUG----------------------   */

   fprintf(stderr, "\n\nPass %d\n", debug++);

   for (i=dstart; i<dend; i++) {  

      fprintf(stderr, "\n%3d: ", i);

      for (j=0; j<3; j++) { fprintf(stderr, " %5ld ", HP[j][i]); }

   }

/*       -----------------------DEBUG -------------------------*/

}



/* 3rd pass  */



i = 2; 

while (i < HPEnd) {

   if ((HP[2][i-1] == 0) && (HP[2][i] - HP[2][i-2] < BasesPerLabel)) { /*  then there was mistake */

      /* and the empty el should then carry the corrected value */

      HP[2][i-1] = HP[2][i-2] + BasesPerLabel;

      HP[0][i-1] = HP[0][i];

      HP[1][i-1] = HP[1][i];

      HP[2][i] = 0; /* and set the bad slot to 0s so we don't get dups */

      HP[0][i] = 0;

      HP[1][i] = 0;

   }

   i++;

}





if (debug > 0) {

/*  --------------------  DEBUG----------------------   */

   fprintf(stderr, "\n\nPass %d\n", debug++);

   for (i=dstart; i<dend; i++) {  

      fprintf(stderr, "\n%3d: ", i);

      for (j=0; j<3; j++) { fprintf(stderr, " %5ld ", HP[j][i]); }

   }

/*       -----------------------DEBUG -------------------------*/

}



  

   /* now go thru all the sites in HP and mark the labels where they're supposed to go.  for each, have to 

      calculate the degree of offset, draw a tic, draw a line to the label, write the label in a small font. */

   fprintf(PS, "/Helvetica findfont 6 scalefont setfont\n");

   NH.radius = PtRadius; 

   NH.xc = PtCenterX; 

   NH.yc = PtCenterY;   /* init NH struct */

   NH.xet = NH.yet = 0;

   NH.alpha = NH.beta =  NH.delta = 0;

   topx = NH.x1 = PtCenterX; 

   topy = NH.y1 = PtCenterY + PtRadius; 





/* ------------------------------ Major and minor tics ---------------------------------------- */

   /* make the major tics -  this should be mechanized, no */

   if (seq_len <= 10000000) {

      if (seq_len <= 10000000 && seq_len > 1000000) { 

   	   tacg_Draw_Seq_Tics(10, 0.5, seq_len, &NH, 1, 1000000,  -1,   "/Helvetica findfont 6 scalefont setfont", PS);

   	   tacg_Draw_Seq_Tics(5,  0.2, seq_len, &NH, 1, 100000,  1000000, "/Helvetica findfont 4 scalefont setfont", PS);

  		   tacg_Draw_Seq_Tics(2,  0.2, seq_len, &NH, 0, 10000,  100000, "/Helvetica findfont 4 scalefont setfont", PS); 

	   } else if (seq_len <= 1000000 && seq_len > 100000) { 

   	   tacg_Draw_Seq_Tics(10, 0.5, seq_len, &NH, 1, 100000,  -1,   "/Helvetica findfont 6 scalefont setfont", PS);

   	   tacg_Draw_Seq_Tics(5,  0.2, seq_len, &NH, 1, 10000,  100000, "/Helvetica findfont 4 scalefont setfont", PS);

  		   tacg_Draw_Seq_Tics(2,  0.2, seq_len, &NH, 0, 1000,  10000, "/Helvetica findfont 4 scalefont setfont", PS); 

      } else if (seq_len <= 100000 && seq_len > 10000) {

   	   tacg_Draw_Seq_Tics(10, 0.5, seq_len, &NH, 1, 10000,  -1,   "/Helvetica findfont 6 scalefont setfont", PS);

   	   tacg_Draw_Seq_Tics(5,  0.2, seq_len, &NH, 1, 1000,  10000, "/Helvetica findfont 4 scalefont setfont", PS);

  		   tacg_Draw_Seq_Tics(2,  0.2, seq_len, &NH, 0, 100,  1000, "/Helvetica findfont 4 scalefont setfont", PS);

      } else if (seq_len <= 10000 && seq_len > 1000) {

   	   tacg_Draw_Seq_Tics(10, 0.5, seq_len, &NH, 1, 1000,  -1,   "/Helvetica findfont 6 scalefont setfont", PS);

		   tacg_Draw_Seq_Tics(5,  0.2, seq_len, &NH, 1, 100,  1000, "/Helvetica findfont 4 scalefont setfont", PS);

  		   tacg_Draw_Seq_Tics(2,  0.2, seq_len, &NH, 0, 10,  100, "/Helvetica findfont 4 scalefont setfont", PS);



      } else {

   	   tacg_Draw_Seq_Tics(10, 0.5, seq_len, &NH, 1, 100,  -1,   "/Helvetica findfont 6 scalefont setfont", PS);

		   tacg_Draw_Seq_Tics(5,  0.2, seq_len, &NH, 1, 10,  100, "/Helvetica findfont 4 scalefont setfont", PS);

      }

   } else {

   	fprintf(stderr, "\nSequence is > 10,000,000; tics jammed too close for this now\n\n");

   }

   

   /* & minor tics */

   i = 1;

/* 

   if (seq_len > 20000) {i = 0;} 

   tacg_Draw_Seq_Tics(5,  0.2, seq_len, &NH, i, 100,  1000, "/Helvetica findfont 4 scalefont setfont", PS);

  	

   if (seq_len < 6000) {

   	tacg_Draw_Seq_Tics(2,  0.1, seq_len, &NH, 0, 10,  100, "/Helvetica findfont 4 scalefont setfont", PS);

   }

 */   

   NH.x1 = PtCenterX; 

   NH.y1 = PtCenterY + PtRadius; 







/* ------------------------------     RE labeling     ---------------------------------------- */



   HPi = 0;

   while (HPi < HPMax) { 

      NH.radius = PtRadius+3; /* '+3' to allow a band to show degeneracies */

      while (HP[0][HPi] == 0 && HPi < HPMax) { HPi++; } /* walk up to the next populated slot */

      if (HPi < HPMax) { /*  we're not to the end, so finish ... */

         RWC = HP[0][HPi];  /* get RWC */

         n = HP[1][HPi];  /* get REi */

   /*       fprintf(stderr, "\nDEBUG:HPi=%ld RWC = %ld RE=%s RWC adj=%ld\n", HPi, RWC, RE[n].E_nam, HP[2][HPi]);  */

         /* calc the coordinates for the next hit position  */ 



         fprintf(PS, "\n\ngsave\n   newpath  "); 

         NH.alpha =  ((double)RWC / (double)seq_len) * 360; /* calc the offset in degrees */



         Calc_Next_Hit_Coords (&NH); /* now the new coords are in NH.xn, NH.yn */

         /* place a ticmark there */

         fprintf(PS, "   %6.2f %6.2f moveto      %% RE#: %d hit=%ld, angle=%f\n", 

                     NH.xn, NH.yn, n, RE[n].Sites[i], NH.alpha);   /* do a direct moveto */

         /* write the tic */

         /* NH.alpha = 0; */

         NH.radius += 10;

         Calc_Next_Hit_Coords (&NH); /* now the new coords are in NH.xn, NH.yn */



         NH.xet = NH.xn;

         NH.yet = NH.yn;

         fprintf(PS, "   1 setlinewidth %6.2f %6.2f lineto stroke   %% draw tic\n", NH.xet, NH.yet);

         if (HP[0][HPi] != 0) {    /* HP[2][HPi]  or put in a very small nonzero */

            NH.alpha =  ((double)(HP[2][HPi]) / (double)seq_len) * 360; /* calc the Label offset in degrees */

            NH.radius += 15;

            Calc_Next_Hit_Coords (&NH); /* now the new coords are in NH.xn, NH.yn */

            /* offset them by the correct amount */

            fprintf(PS, "   %6.2f %6.2f moveto 0.1 setlinewidth  %6.2f %6.2f lineto stroke  %% angle=%f\n", 

                           NH.xet, NH.yet, NH.xn, NH.yn, NH.alpha );         

            fprintf(PS, "   %6.2f %6.2f  moveto   %6.2f rotate  ", 

            NH.xn, NH.yn, (90 - NH.alpha));

         }   



         /* compose & print the label */

         if (RE[n].E_olap > 0)      { o = "\\\\";}

         else if (RE[n].E_olap > 0) { o = "/";}

         else                       { o = "|";}

         star = ' ';

         if (RE[n].E_Ncuts == 1)    { star = '*';}

         sprintf(Label, "%s [%ld] %s %c", RE[n].E_nam, RWC, o, star);



         if (NH.alpha > 180) {         

            fprintf(PS, " gsave\n 180 rotate\n ");                                         /*  */

            fprintf(PS, " (%s) stringwidth neg exch neg exch rmoveto\n (%s) show\n grestore \n grestore\n \n" , Label, Label);

         } else { 

            fprintf(PS, "(%s) show\ngrestore\n", Label); /*  show the label in correct position */

         }

         HPi++;

      }

   } 

   

   

   /* ------------------------------       Plot ORFs       ---------------------------------------- */

   

   /* now if ORFs were requested, print out the ORFS greater than the requested minimum */

   fprintf(PS, "\n gsave\n");

   NH.radius = PtRadius - 30; /* start point for the ORF arcs */

   ORFNum = 0;

   ORF_density = 0.0;

   FCnt = 0;

   ORF_line_width = 3; /* temp est  */

   

   /* this has to be modified to do 1-3 1st and 3-6 after, compacting output so that if 1246, it will draw

      1,2 then go directly to 4 then to 6, eliminating circles for 3,5. Do this by sorting the ORFs and compact 

      to be contiguous.  Also have to reverse the arrowhead symbol for frames 3-6  */

      

   fprintf(PS, "\n\n/Helvetica findfont 8 scalefont setfont\n");

   while (F.ORFs2Do[FCnt] != 0) { /* do all at first, then follow the F.ORFs2Do array as to which to print */

      f = F.ORFs2Do[FCnt] - 1;

      ORFNum = 0;

      NH.radius -= 10 ; /* make radius a bit smaller for each frame. */

      hairline_radius = NH.radius + 7;

      /* draw the hairline circle as a separator */

      NH.alpha = 0;

      Calc_Next_Hit_Coords (&NH); 

      if (f<3) { arrow = ">"; }

      else     { arrow = "<"; }

      fprintf(PS, "\n newpath %6.2f %6.2f moveto  (%d%s) show \n", NH.xn, (NH.yn-2), (f+1), arrow); 

      fprintf(PS, "\n newpath 0.1 setlinewidth 0.5 setgray\n %f %f %f  0 360  arc stroke\n", 

                    NH.xc, NH.yc, hairline_radius-2 );

      ORF_density += 0.07;

      /* ORF_struct ORFs contains 6 elements, populated in order..? that contain the ORFs found.

         at each ORF iter, have to move to the appro radius and draw the arcs with the appro shading

         or hatching (identifying table at the side*/

      while (ORFs[f][ORFNum].orflength > 0) { /* indicator that there are no more ORFs in this frame */

         /* calc ORF angular offsets in degrees for the beginning of the ORF */

         if (f > 2) {

            temp = ORFs[f][ORFNum].E_offsetBP;

            ORFs[f][ORFNum].E_offsetBP = ORFs[f][ORFNum].B_offsetBP;

            ORFs[f][ORFNum].B_offsetBP = temp;

         }

         NH.alpha =  ((double)ORFs[f][ORFNum].E_offsetBP / (double)seq_len) * 360; /* calc the offset in degrees */

         Calc_Next_Hit_Coords (&NH); /* now the new coords are in NH.xn, NH.yn */     

         arc_end = NH.alpha;    

         /* and move there */

         fprintf(PS, "\n %f %f moveto\n", NH.xn, NH.yn); /* unless we have to reset to the origin 1st  */

         /* possibly add an arrowhead here to show direction of transcription/translation */

         NH.alpha =  ((double)ORFs[f][ORFNum].B_offsetBP / (double)seq_len) * 360; /* same thing with the End */

         Calc_Next_Hit_Coords (&NH); /* now the new coords are in NH.xn, NH.yn */  

         /* now write the arc from the ORF origin to the end at line width and density set already */       

         fprintf(PS, "%d setlinewidth %f setgray\n newpath  %f %f %f  %f %f  arc "

         "%% frame=%d ORF=%d begin=%ld end=%ld len=%d\n  stroke\n", 

                     ORF_line_width, ORF_density,  NH.xc, NH.yc, NH.radius, 90-arc_end, 90-NH.alpha, 

            f, ORFNum, ORFs[f][ORFNum].B_offsetBP, ORFs[f][ORFNum].E_offsetBP, ORFs[f][ORFNum].orflength);

         

         ORFNum++;

      } 

      FCnt++;  

   }

   /* and plot a final hairline to finish off the ORF */

   NH.radius -= 4 ; /* make radius a bit smaller  */

   NH.alpha = 0;

   Calc_Next_Hit_Coords (&NH); 

   fprintf(PS, " newpath 0 setgray 0.2 setlinewidth\n %f %f %f  0 360  arc stroke\n", 

                    NH.xc, NH.yc, NH.radius);

    

/* ------------------------------     Plot degeneracies    ---------------------------------------- */

	

   fprintf(PS, "\n\n%% Marking degeneracies\n3 setlinewidth 0.0 setgray\n");

   if (F.LogDegens) {

      NH.radius = PtRadius +1.5;

      NH.x1 = NH.xc;

      NH.y1 = NH.yc + NH.radius;

      for (i=2; (i<Degen_Log[1] && Degen_Log[i] < seq_len); i+=2) {

         NH.alpha = ((double)Degen_Log[i] / (double)seq_len) * 360; /* calc the offset in degrees to the START of the degen run */

         

         Calc_Next_Hit_Coords (&NH);

         arc_end = NH.alpha;  

         fprintf(PS, "\n newpath %6.2f %6.2f moveto newpath %% Base %ld \n", 

                               NH.xn, NH.yn,                   Degen_Log[i]);

         NH.alpha = ((double)Degen_Log[i+1] / (double)seq_len) * 360; /* offset to END of the degen run */

         fprintf(PS, "  %f %f %f  %f %f  arc stroke\n\n", NH.xc, NH.yc, NH.radius, 90-NH.alpha, 90-arc_end ); /* 90-NH.alpha, 90-arc_end */

      }

   }

   

   fprintf(PS, "\nshowpage\n");

   fclose(PS);

}



   /* int tic_size - the length of the tics in pts

      float line_width - the width of the line in pts

      long seq_len - the sequence length in base pairs

      struct NextHit_struct *NH - the Next Hit struct 

      int drawlabels - whether to print labels or just write the tics

      long repeat_period - the tic period (1000 -> 1000, 2000, 3000, etc)

      int skip_period - what period should be skipped.  If already done repeat_period -> 1000 and

            doing minor tics at 100, don't want to overwrite the 1000 label 

      int font_string - the font characteristics as Postscript likes it:

         

      FILE *PS - the file pointer to the file we're writing to

   */    

void tacg_Draw_Seq_Tics(int tic_size, float line_width, long seq_len, struct NextHit_struct *NH, int drawlabels,

                        long repeat_period, int skip_period, char *font_string, FILE *PS ){

   int i, ntics, nk; 

/*   label_offset = -35;  how far should we move towards the center to compensate for the label */

   char *Label;

   if ((Label = calloc(128,sizeof(char))) == NULL) { BadMem("tacg_Draw_Seq_Tics:Label", 1); }



   tic_size *= -1;  /* tics go towards center */

   ntics = (int)(seq_len / repeat_period);

/*   if ((seq_len % 1000) < 0.00001) { ntics--; }  knock off the last tic if it overlaps with the origin. */

   nk=0; /* num of k's, 1000s were the orig tic counter */

   fprintf(PS, "gsave\n %f setlinewidth  %s\n", line_width, font_string);

   for (i=0; i<= ntics; i++) {  /*  && ((*NH).alpha < 359) */

      (*NH).alpha = (double)nk / (double)seq_len * 360;

      if ((*NH).alpha < 360 ) { /* && (*NH).alpha > 0 */

      fprintf(PS, "\n\nnewpath\n gsave\n ");  

      Calc_Next_Hit_Coords (NH);

      

      fprintf(PS, " %f %f moveto \t%% tic mark: %d, angle=%f\n", 

                        (*NH).xn, (*NH).yn, nk, (*NH).alpha); 

      fprintf(PS, " %f rotate\n gsave\n", (90-(*NH).alpha));        /* rotate */

      

      if (nk > 0) sprintf(Label, "%d", nk);

      else if (skip_period < 0 ) sprintf(Label, "     0 "); 

      else sprintf(Label, " ");  

      

      fprintf(PS, " %d 0 rlineto \n stroke\n grestore\n", tic_size);



      if (drawlabels) {  /* for the TICS */

         if ((*NH).alpha > 180) {         

            fprintf(PS, " gsave\n 180 rotate\n ");

            fprintf(PS, " (%s) stringwidth neg exch .5 mul exch  rmoveto\n (%s) show\n grestore\n grestore\n" , Label, Label);

         } else { 

            fprintf(PS, " (%s) stringwidth neg exch dup 2.5 mul sub exch  rmoveto\n",  Label);

            fprintf(PS, "(%s) show\n grestore\n\n", Label); /*  show the label in correct position */

         }

      } else { fprintf(PS, " grestore\n"); }    

      nk += repeat_period;

      if (skip_period > 0 && nk % abs(skip_period) < 0.00001) { nk += repeat_period; } /* skip the skip period */

   } 

   }

    fprintf(PS, " grestore\n");    /* to restore the context for whatever was set before */

} 





/* passed the pointer of the populated struct, it calculates the x,y coords of the next hit and passes them

   back as struct elements xn,yn */

void Calc_Next_Hit_Coords (struct NextHit_struct *NH) {

   double pi = 3.141592653589, reflect = 1;

   if ((*NH).alpha > 180) {reflect = -1; }

   (*NH).x1 = (*NH).x1 - (*NH).xc;

   if ( fabs((*NH).x1) > (*NH).radius) {

/*       fprintf(stderr, "abs(x1) (%f) > radius (%f)?? \n", (*NH).x1, (*NH).radius);  */

      (*NH).x1 = (*NH).radius; /* make sure it gets reset before next iteration */

   } 

   (*NH).y1 = (*NH).y1 - (*NH).yc;

   (*NH).beta =  (pi * asin((*NH).x1 / (*NH).radius)) / 180; 

   (*NH).delta = 90 - (*NH).alpha - (*NH).beta;

/*    fprintf(stderr, "\n\tbeta=%f delta=%f alpha=%f asin=%f\n", 

      (*NH).beta, (*NH).delta, (90 - (*NH).beta - (*NH).delta), asin((*NH).x1 / (*NH).radius)); */

   (*NH).yn = (sin((*NH).delta / (360 / (2 * pi))) * (*NH).radius); /* don't re-add the center value yet */

   (*NH).xn = (*NH).xc + (sqrt((*NH).radius * (*NH).radius - (*NH).yn * (*NH).yn) * reflect);

   (*NH).yn = (*NH).yc + (*NH).yn; /* now add it back */

   

   (*NH).x1 = (*NH).xn; (*NH).y1 = (*NH).yn;    /* update the 'old' coords as well for the next round */

}