#!/usr/bin/perl

# read in 4 xpk files (dif and sum of aligned and unaligned sample)
# constructed from IPAP data (2D spectra and xpk files)
# and calculate RDC for assigned peaks

# calculates aligned minus isotropic.
# if one of the atoms is N then the sign is reversed.
# so the sign of the RDC should be correct.

# option for normalizing DC to NH vector

# CAUTION - the sign of hnco DC may need to be obtained elsewhere.
# Do not rely on this script for hnco.

# RDC scaling factors relative to N-H are
#    N-H      1.0
#    C-H     -2.1
#    Co-Ca   -1/5.05 = -0.198
#    HN-Co   -1/3.2  = -0.312
#    N-Co     1/8.33 =  0.120

# added type #5 for NhsqcS3 sequence 2007

#print "unaligned diff xpk: ";
#chop($u_dif = <STDIN>);
#print "unaligned sum xpk: ";
#chop($u_sum = <STDIN>);
#print "aligned diff xpk: ";
#chop($a_dif = <STDIN>);
#print "aligned sum xpk: ";
#chop($a_sum = <STDIN>);
$u_dif="udiflist.xpk";
$u_sum="usumlist.xpk";
$a_dif="adiflist.xpk";
$a_sum="asumlist.xpk";

print "1. N-Hn ipap (2D or 3D)\n";
print "2. N-Co from hnco\n";
print "3. Hn-Co from hnco\n";
print "4. Co-Ca from hnco or 2D-S3\n";
print "5. N-Hn from NhsqcS3 (coupling in H dimension)\n";
print "Select type of experiment: ";
chop($type = <STDIN>);
if($type == 2) {
    print "enter value of Jf: ";
    chop($Jf = <STDIN>);
    }
if($type == 4) {
    print "enter lambda scaling for 2D-S3, if any: ";
    chop($lambda = <STDIN>);
    if($lambda == 0) {
        $lambda=1; }
    }
$norm = 0;
if ($type >= 2) {
	print "normalize to NH ?: ";
	chop($_ = <STDIN>);
	if(/y/i) {
		$norm=1; }
	}

print "output RDC file name: ";
chop($output = <STDIN>);
$udi=0;
$usi=0;
$adi=0;
$asi=0;
open(IN,$u_dif) || die "open unaligned diff failed";
for ($i=0; $i<4; $i++) {
    <IN>;
    }
$_=<IN>;
s/\{//g;
s/\}//g;
@w=split;
if($type == 3 || $type==5) {
    $freq = $w[0];
    } else {
            $freq = $w[1];
            }
print "frequency is $freq \n";
$_=<IN>;
while (<IN>) {
    chop;
    @w = split;
    if ($w[1] ne "{?}" && $w[8] ne "{?}") {
        $ud1[$udi] = $w[1];
        $ud2[$udi] = $w[8];
        $uds1[$udi] = $w[2];
        $uds2[$udi] = $w[9];
        $udi++;
        }
    }
close(IN);
open(IN,$u_sum) || die "open unaligned sum failed";
for ($i=0; $i<6; $i++) {
    <IN>;
    }
while (<IN>) {
    chop;
    @w = split;
    if ($w[1] ne "{?}" && $w[8] ne "{?}") {
        $us1[$usi] = $w[1];
        $us2[$usi] = $w[8];
        $uss1[$usi] = $w[2];
        $uss2[$usi] = $w[9];
        $usi++;
        }
    }
close(IN);
open(IN,$a_dif) || die "open aligned diff failed";
for ($i=0; $i<6; $i++) {
    <IN>;
    }
while (<IN>) {
    chop;
    @w = split;
    if ($w[1] ne "{?}" && $w[8] ne "{?}") {
        $ad1[$adi] = $w[1];
        $ad2[$adi] = $w[8];
        $ads1[$adi] = $w[2];
        $ads2[$adi] = $w[9];
        $adi++;
        }
    }
close(IN);
open(IN,$a_sum) || die "open aligned sum failed";
for ($i=0; $i<6; $i++) {
    <IN>;
    }
while (<IN>) {
    chop;
    @w = split;
    if ($w[1] ne "{?}" && $w[8] ne "{?}") {
        $as1[$asi] = $w[1];
        $as2[$asi] = $w[8];
        $ass1[$asi] = $w[2];
        $ass2[$asi] = $w[9];
        $asi++;
        }
    }
close(IN);
#for ($i=0; $i<$udi; $i++) {
#    print "$i   $ud1[$i]    $ud2[$i] \n";
#    }

$u=0;
for ($i=0; $i<$udi; $i++) {
    for ($j=0; $j<$usi; $j++) {
        if ($ud1[$i] eq $us1[$j] && $ud2[$i] eq $us2[$j]) {
            if($type == 3) {
                $unalign[$u] = abs($uds1[$i] - $uss1[$j]);
            } else {
                $unalign[$u] = abs($uds2[$i] - $uss2[$j]);
                }
            $unalignN[$u] = $us2[$j];
            $unalignH[$u] = $us1[$j];
#print "$u  $unalignN[$u]   $unalign[$u]    $uds2[$i] - $uss2[$j] \n";
            $u++;
            }
        }
    }
$a=0;
for ($i=0; $i<$adi; $i++) {
    for ($j=0; $j<$asi; $j++) {
        if ($ad1[$i] eq $as1[$j] && $ad2[$i] eq $as2[$j]) {
            if($type == 3) {
                $align[$a] = abs($ads1[$i] - $ass1[$j]);
            } else {
                $align[$a] = abs($ads2[$i] - $ass2[$j]);
                }
            $alignN[$a] = $as2[$j];
            $alignH[$a] = $as1[$j];
#print "$a  $alignN[$a]   $align[$a]    $ads2[$i] - $ass2[$j] \n";
            $a++;
            }
        }
    }

$err = 0.5;
$order = 1.0;
$count=0;
open(OUT,">$output");
for ($i=0; $i<$u; $i++) {
    for ($j=0; $j<$a; $j++) {
        if ($unalignN[$i] eq $alignN[$j]) {
            $rdc[$count] = ($align[$j] - $unalign[$i]) * $freq;
            $temp = index($alignN[$j],".") - 1;
            $temp2 = index($alignN[$j],"}") - 2;
            $rdcN1[$count] = substr($alignN[$j],1,$temp);
            $rdcA1[$count] = substr($alignN[$j],$temp+2,$temp2-$temp);
            $temp = index($alignH[$j],".") - 1;
            $temp2 = index($alignH[$j],"}") - 2;
            $rdcN2[$count] = substr($alignH[$j],1,$temp);
            $rdcA2[$count] = substr($alignH[$j],$temp+2,$temp2-$temp);
            if($type == 2) {
                $rdcA2[$count] = "C";
                $rdcN2[$count] = $rdcN1[$count] - 1;
                $rdc[$count] = $rdc[$count] / ($Jf+1);
			if($norm == 1) {
				$rdc[$count] = $rdc[$count] * 8.33; }
			}
            if($type == 3) {
                $rdcA1[$count] = "C";
                $rdcN1[$count] = $rdcN2[$count] - 1;
			if($norm == 1) {
				$rdc[$count] = $rdc[$count] * (-3.2); }
                }
		if($type == 4) {
			$rdcA2[$count] = "CA";
                if($rdcA1[$count]="N") {    # fix for 2D N-NH of S3 exp
			    $rdcA1[$count] = "C";
                    $rdcN1[$count] = $rdcN1[$count]-1;
                    }
			$rdcN2[$count] = $rdcN1[$count];
                $rdc[$count] = $rdc[$count] / $lambda;
			if($norm == 1) {
				$rdc[$count] = $rdc[$count] * (-5.05); }
			}
            if ($type == 1 || $type == 2 || $type ==5) {
                $rdc[$count] = -$rdc[$count]; }
            printf OUT "%-3i %3s %3i %3s  %8.3f    %3.1f    %4.2f\n",$rdcN1[$count], $rdcA1[$count], $rdcN2[$count], $rdcA2[$count], $rdc[$count], $err, $order;
            $count++;
            }
        }
    }

print "$u unaligned and $a aligned crosspeaks \n";
print "$count matching peaks found \n";