#!/Users/gunnar/bin/yorick -batch

/* CONVERT ksdssp to direx format */
/* USAGE:  ksdssp2dx file.pdb ksdssp.out */

#include "string.i"


main
{

args=get_argv();

if (numberof(args()) != 3) {
 print,"USAGE: ksdssp2dx file.pdb ksdssp.out "
 exit;
}

fname=args(2);
f2name=args(3);

f=open(fname) 
old_resid=-1;
n=0;
rid_number=-1;
while ( line=rdline(f) ) {
 first_col=strpart(line,1:4);
 if ( first_col=="ATOM"  || first_col=="HETATM" ) {
     rid=strpart(line,23:26);
     sread,rid,format="%i",rid_number
     if (old_resid != rid_number)  {
         old_resid= rid_number;
         n++ ;
     }
 }
}
close,f;
print,"number of residues",n;


chain=array("A",n);
resid=array(int,n);
sec=array(" ",n);

f=open(fname);
n=0;
old_resid=-1;
rid_number=-1;
while ( line=rdline(f) ) {
  if ( strpart(line,1:4)=="ATOM" ) {
    rid=strpart(line,23:26);
    sread,rid,format="%i",rid_number
    if (old_resid != rid_number) {
       resid(++n)=rid_number;     
       chain(n)=strpart(line,22:22);
       old_resid= rid_number;
    }
  }
}
close,f;

istart=-1;
istop=-1;
f=open(f2name)
while ( line=rdline(f) ) {
  if ( strpart(line,1:5)=="HELIX" ) {
   istart_str=strpart(line,21:25);
   istop_str=strpart(line,33:37);
   sread,istart_str,format="%i",istart
   sread,istop_str,format="%i",istop
   chain_start=strpart(line,20:20)
   chain_stop=strpart(line,32:32)
   if (chain_start == chain_stop) {
     print,chain_start
     for (i=istart;i<=istop;i++)   {
       ind=where(chain==chain_start)
       ind2=where(resid(ind)==i)(0)
       //print,i,resid(ind(ind2))
       //print,i,ind(ind2),resid(ind(ind2)),chain(ind(ind2))
       //print,i
       sec(ind(ind2))="H";
     }
   }
  }

  if ( strpart(line,1:5)=="SHEET" ) {
    istart_str=strpart(line,23:26);
    istop_str=strpart(line,34:37);
    sread,istart_str,format="%i",istart
    sread,istop_str,format="%i",istop
    chain_start=strpart(line,22:22)
    chain_stop=strpart(line,33:33)
    if (chain_start == chain_stop) {
     for (i=istart;i<=istop;i++)   {
       ind=where(chain==chain_start)
       ind2=where(resid(ind)==i)(0)
       sec(ind(ind2))="E";
     }
    } 
  }

}
close,f



print, "really done"

fout=open("out_secstr.dat","w")
for (i=1;i<=n;i++) {
  write(fout,format="%i  %i   %s  X  %s  \n",i,resid(i),chain(i),sec(i))
}
close,fout
 

}