#!/usr/bin/perl

# calculate the change in chemical shifts of N and HN
# from xpk files of 2 Nhsqc spectra
# Nscale is scaling factor for N15 relative to H

$Nscale=0.15;
print "compute difference of Nhsqc chemical shifts\n";
print "xpk #1: ";
chop($u_dif = <STDIN>);
print "xpk #2: ";
chop($a_dif = <STDIN>);


print "output file name: ";
chop($output = <STDIN>);
$udi=0;
$adi=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) {
#    $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,$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);
#for ($i=0; $i<$udi; $i++) {
#    print "$i   $ud1[$i]    $ud2[$i] \n";
#    }


$count=0;
open(OUT,">$output");
for ($i=0; $i<$udi; $i++) {
    for ($j=0; $j<$adi; $j++) {
        if ($ad1[$j] eq $ud1[$i] && $ad2[$j] eq $ud2[$i]) {
            $temp = index($ud1[$i],".") - 1;
            $temp2 = index($ud1[$i],"}") - 2;
            $resn[$count] = substr($ud1[$i],1,$temp);
            $atom[$count] = substr($ud1[$i],$temp+2,$temp2-$temp);
            $dif1 = $ads1[$j] - $uds1[$i];
            $dif2 = $ads2[$j] - $uds2[$i];
            $diff[$count] = sqrt($dif1*$dif1 + $Nscale*$Nscale*$dif2*$dif2);
            printf OUT "%-3i %6.3f\n", $resn[$count],$diff[$count];
            $count++;
            }
        }
    }

print "$udi reference and $adi sample crosspeaks \n";
print "$count matching peaks found \n";