#!/usr/bin/awk -f
#Notes:
#We may want to imlement a tolerance for small deviations.
#We might want to create a second list of parameters (hits),
#sort them and print them. Parts of the code are shaped to
#accomodate this functionality. However, we'll have to re-introduce
#the qsort and qlookup functions that were (partially) present
#in version 0.7. The present version branched off to save time before
#any debugging took place (ie. 0.7 is full of bugs).
BEGIN {
#initialize constants
  j=32;
  priorities="CGNGOGPGSGFGCLBRIGHELIBEBGNENAMGALSIARKGCASCTIVGCRMNFECONICUZNGAGEASSEKRRBSRYGZRNBMOTCRURHPDAGCDINSNSBTEXEHGX ";
  for(i=1;i<110;i+=2) prior[substr(priorities,i,2)]=sprintf("%c",j++);
  c255=sprintf("%c",255);
  hel="";
  for(i=1;i<=22;i++) hel=hel c255;
#read parameter files
  parf[1]=par1;parf[2]=par2;
  for(i in parf) {
    parfile=parf[i];npp="";recnum=0;
    do {gl=getline < parfile;
      recnum++;
      if ($0 ~ "^ANGLES") npp=4;
      if ($0 ~ "^DIHEDRALS") break;
      if (npp) {
        if ($1!="" && $4!="" && $1!~"^!") {
	  rec[1]=substr($0,1,6);
	  rec[2]=substr($0,8,6);
	  rec[3]=substr($0,15,6);
          rec[4]=$5;
	  rec[5]=parfile;
	  rec[6]=recnum;
	  parpri(1) } } }
     while (gl) } ;
#initialize values for main program
  rsel=toupper(rsel)
  if (fp) {
    rec[1]=substr(toupper(fp1) "     ",1,6);
    rec[2]=substr(toupper(fp2) "     ",1,6);
    rec[3]=substr(toupper(fp3) "     ",1,6);
    parpri(0);
    tofind=ident }
  first=1}

{$0=toupper($0);
if ($1=="RESI" || $1=="PRES" || $1=="END") {
  if (first) first=""
   else if (reslin~rsel) {
#core functionality
    ntri=0;
  #detect trigonal centers. The actual summation is done in a seperate loop
  #to allow future extension.
    for (i in con) if (con[i]==2) trig[++ntri]=i;
  #sum over rings
    for (i=1;i<=nrings;i++) {
      sumang=0;
      found="";
      j2=0;j3=1;max=rmem[i];
      for(j=0;j<max;j++) {
	matchpar(ratom[i,j],ratom[i,(++j2)%max],ratom[i,(++j3)%max]);
#print j " " (++j2) " " (++j3)%max " --- " ratom[i,j] " " ratom[i,j2%max] " " ratom[i,j3%max];
#matchpar(ratom[i,j],ratom[i,j2%max],ratom[i,j3%max]);
#print ident;
        inv[j]=ident};
      tarang=(max-2)*180;
      if (!fp && (sumang < tarang-tol || sumang > tarang+tol)) found=1;
      if (found) {
	printf ("%s ring %s : sum = %s (%i)\n",res,rspec[i],sumang,tarang);
	verbose(max) } };
  #sum over trigonal centers
    for (i=1;i<=ntri;i++) {
      sumang=0;
      found="";
      c=trig[i];
      b1=bat[c,0];b2=bat[c,1];b3=bat[c,2];
      matchpar(b1,c,b2);inv[0]=ident;
      matchpar(b2,c,b3);inv[1]=ident;
      matchpar(b3,c,b1);inv[2]=ident;
      if (!fp && (sumang < 360-tol || sumang > 360+tol)) found=1;
      if (found) {
	printf ("%s center %s bound to %s %s %s : sum = %s (360)\n",
res,c,b1,b2,b3,sumang);
        verbose(3) } } };
#initialize residue data
  res=$1 " " $2;
  reslin=substr($0,5);
  delete atype; delete con;
  nrings=0 } ;
#read in residue data
if ($1=="!RING") {
  rmem[++nrings]=$2;
  rspec[nrings]=substr($0,7);
  arg=2;
  for(i=0;i<$2;i++) ratom[nrings,i]=$++arg }
if ($1=="ATOM") {atype[$2]=substr($3 "     ",1,6);con[$2]=-1}
if ($1=="BOND" || $1~"^DOUB") {
  arg=1;
  while ($++arg!="") {
    at1=$arg; at2=$++arg;
    c1=++con[at1];c2=++con[at2];
    bat[at1,c1]=at2;bat[at2,c2]=at1 } } }

function verbose(mm) {
  print "Parameters involved:";
  for(n=0;n<mm;n++) {
    id=inv[n];
    printf ("%s line %i: %s %s %s ang=%s\n",
par[5,id],par[6,id],par[1,id],par[2,id],par[3,id],par[4,id]) } }

function matchpar(an1,an2,an3) {
  rec[1]=atype[an1];rec[2]=atype[an2];rec[3]=atype[an3];
  parpri(0);
  if (fp && ident==tofind) {
    printf("Parameter found in %s: %s %s %s\n",res,an1,an2,an3);
    found=1 }
  sumang+=par[4,ident] }

#We could move the "add" functionality to the main program but choose to keep
#it here to allow future extension.
#Also, it would be good practise to make parpri return ident
#(and zero in case of nocom)
function parpri(add) {
  nocom=1;
  for(n=1;n<=3;n++) {
    ri=rec[n];
    pri=prior[substr(ri,1,2)];
    if (pri=="") {nocom="";break};
    att[n]=pri substr(ri,3)}
  if (nocom) {
    if (att[1]>att[3]) swap(1,3);
    ident=att[2] att[1] att[3];
    if (add) for(n=1;n<=6;n++) par[n,ident]=rec[n] } }

function swap(x,y) {
  z=att[x];att[x]=att[y];att[y]=z;
  z=rec[x];rec[x]=rec[y];rec[y]=z }
