#!/usr/bin/perl

# convert RDC file into format for pales and ehm
# for ehm use normalized RDC data

print "input rdc file name: ";
$fname = <STDIN>;
open(INFILE, $fname);
print "input pdb file name(leave blank for .seq): ";
$fname2 = <STDIN>;
chop($fname2);
$ispdb = 1;
if ($fname2 eq "") {
	print "input sequence file name: ";
	$fname2 = <STDIN>; 
	$ispdb = 0; }
open(INFILE2, $fname2);
print "output pales rdc file name: ";
$fname3 = <STDIN>;
$n=0;
while (<INFILE>) {
        chop;
        ($resn1[$n],$atom1[$n],$resn2[$n],$atom2[$n],$rdc[$n],$err[$n],$wt[$n]) = split;
#       $rdc[$num] = -$rdc[$num];
        $n++;
    }
close(fname);
#for ($i=0; $i<$n; $i++) {
#        $resn[$n+$i] = $resn[$i]+100;
#        $rdc[$n+$i] = $rdc[$i];
#        }
#$n=$n*2;

$i=0;
while (<INFILE2>) {
	if ($ispdb == 1) {
		$key = substr($_, 0,4);
        if ($key eq "ATOM") {
        # dissect line into parts
        chop;
        $num = substr($_,6,5);
        $atom = substr($_,12,4);
        $res = substr($_,17,3);
        $chainid = substr($_,21,1);
        $resn = substr($_,22,4);
        $xcor = substr($_,30,8);
        $ycor = substr($_,38,8);
        $zcor = substr($_,46,8);
        $occ = substr($_,54,6);
        $bfac = substr($_,60,6);
        $segid = substr($_,72,4);
        $elm = substr($_,76,2);
        $charge = substr($_,78,2);


            if ($i == 0) {
                $resnum[$i] = $resn;
                $resnam[$i] = $res;
#printf "0-> %d %d %d %s %s\n", $i, $resn, $resnum[$i], $res, $resnam[$i];
                $i++;
                }
            else {
                if ($resn > $resnum[$i-1]) {
                    $resnum[$i] = $resn;
                    $resnam[$i] = $res;
#printf "-> %d resn %s resnum %s res %s resnam %s\n", $i, $resn, $resnum[$i], $res, $resnam[$i];
                    $i++;
                    }
                }
            }
	} else {
		chop;
		$resnam[$i] = $_;
		$resnum[$i] = $i+1;
#printf "%s\n", $resn;
		$i++;
		}
	}
close(fname2);
$natoms=$i;

open(OUT,">$fname3");

# print anything here since it usually is not used anyway
print OUT "DATA SEQUENCE AKKKSKLEIIQAILEACKSGSPKTRIMYGANLSYALTGR\n";
print OUT "DATA SEQUENCE YIKMLMDLEIIRQEGKQYMLTKKGEELLEDIRKFNEMRK\n";
print OUT "DATA SEQUENCE NMDQLKEKINSVLSIRQ\n";

print OUT "\nVARS   RESID_I RESNAME_I ATOMNAME_I RESID_J RESNAME_J ATOMNAME_J D DD W\n";
print OUT "FORMAT %5d %6s %6s %5d %6s %6s %9.3f %9.3f %5.2f \n\n";
for ($i=0; $i<$n; $i++) {
	for ($j=0; $j<$natoms; $j++) {
		if ($resn1[$i] == $resnum[$j]) {
			$t1 = $resnam[$j]; }
		if ($resn2[$i] == $resnum[$j]) {
			$t2 = $resnam[$j]; }
		}
#$err[$i]=0.5;
#$wt[$i]=1;
if($atom2[$i] eq "HN") {
	$wt[$i] = 1.0;
	$err[$i] = 1.0; }
if($atom2[$i] eq "C") {
	$wt[$i] = 5.05;
	$err[$i] = 1/5.05; }
if($atom2[$i] eq "CA") {
	$wt[$i] = 8.33;
	$err[$i] = 1/8.33; }
	printf OUT "%5d %6s %6s %5d %6s %6s %9.3f %9.3f %5.2f\n",
	$resn1[$i],$t1,$atom1[$i],$resn2[$i],$t2,$atom2[$i],$rdc[$i],$err[$i],$wt[$i];
	}
close($fname3);