#!/usr/bin/perl

#perl script to get the ratio of diagonal to crosspeak
#intensities/volumes

#can use volumes or intensities, but they must be stored in the
#volume entity in felix, and export to a vol file. There is a Felix
#macro (relvolume.mac) that is supposed to be able to calculate accurate
#peak heights and store them in the volume entity. However, I could
#not get this to work, so I used volumes

#For more information about the calculation for HNHA data, see:
#Vuister and Bax (1993) JACS 115, 7772-7777

#Melanie Nelson
#4/15/99

use Math::Trig;

if ($#ARGV != 1){
	$prg_name = $0;
	$prg_name = `basename $prg_name`;
	chomp ($prg_name);
	die<<EODIE;
This script reads the xpk and vol files exported from Felix97, and
calculates the ratio of the intensity of the diagonal peak (HN:HN) 
and the crosspeak (HN:HA) for all residues, using equation 2
Vuister and Bax (1993) JACS 115, 7772-7777

xpk format:
1   359.950     1.732    15  1:SER_2:HN          53.302 1.184    15  1:SER_2:HN          49.785     0.587     0  1:SER_2 1.00
2   359.984     1.978    15  1:SER_2:HN         221.063 1.584    12  1:SER_2:HA          49.941     0.891   454  1:SER_2:N 1.00

vol format:
     1     0.510E+09     0.000E+00     0.000E+00     0.000E+00 0.000E+00     0.000E+00
     2    -0.944E+08     0.000E+00     0.000E+00     0.000E+00 0.000E+00     0.000E+00

The volume or intensity information to be used should be stored in
the first slot of the volume entity. The xpk and vol files must have
matching peak numbers.

Usage: $prg_name fileroot outputfile
fileroot = the root of the file name for the xpk and vol files (must
be the same, and the two files must have the .xpk and .vol
extensions, respectively)
EODIE
}

#read in xpks and store relationship between assignment and peak
#number

$xpkfile = "$ARGV[0].xpk";
$volfile = "$ARGV[0].vol";

open (XPKFILE, $xpkfile) ||
	die ("Could not open xpk file, $xpkfile\n");

while (<XPKFILE>){

#skip comment and unnecessary lines
	next if /\#/;
	next if /c\*\*/;
	next if /^xpk/;

	@line = split;

#extract atom and residue info
	($junk, $res1, $atom1) = split(/:/,$line[4]);
	($resname1, $resnum1) = split(/_/,$res1);
	($junk, $res2, $atom2) = split(/:/,$line[8]);
	($resname2, $resnum2) = split(/_/,$res2);

#is it a crosspeak or a diagonal peak?
#store the peak number, for matching with volumes/intensities from
#the vol file
	if ($line[4] eq $line[8]){
		$diagpk{$resnum1} = $line[0];
		$resname{$resnum1} = $resname1;
	} else {
		$crosspk{$resnum1}{$atom2} = $line[0];
	}
			
}
close XPKFILE;

#now read in vol file and get volumes for each peak
open (VOLFILE, $volfile) ||
	die ("Could not open vol file, $volfile\n");

while (<VOLFILE>){

#skip comment and unnecessary lines
	next if /\#/;
	next if /c\*\*/;
	next if /^xpk/;

	@line = split;

	$vol{$line[0]} = $line[1];

}
close VOLFILE;

#find out the additional parameters for the equation:
print "What is the T2 delay (from the pulse program, in seconds)?\n";
chomp($tau = <STDIN>);
print "Do you wish to include the correction factor used by Vuister
and Bax (multiply J values by 1.1)? Enter Y or N:\n";
chomp($mult_test = <STDIN>);

open (OUTFILE, ">$ARGV[1]") ||
	die ("Could not open output file $ARGV[1]\n");

$pi = atan2(1,1) * 4;

print OUTFILE "J-couplings from getratio. tau = $tau sec\n";
if ($mult_test eq "Y" or $mult_test eq "y"){
	print OUTFILE "J-couplings are corrected by a factor of 1.1\n";
} else {
	print OUTFILE "J-couplings are not corrected for systematic error\n";
}
print OUTFILE "----------------------------------------------\n";
print OUTFILE "residue   atoms   J-coupling\n";
print OUTFILE "----------------------------------------------\n";

#sort the keys of the diagpk hash, so that output will be in order
sub numerically { $a <=> $b }
@reslist = sort numerically (keys %diagpk);

#foreach $x (keys %diagpk){
foreach $x (@reslist){
	$dpk = $diagpk{$x};
	$dpk_vol = $vol{$dpk};
	foreach $atm ( keys %{ $crosspk{$x} } ){
		$cpk = $crosspk{$x}{$atm};
		$cpk_vol = $vol{$cpk};

#this step takes into account the difference in sign for crosspeaks
#that have been folded
		if ($cpk_vol > 0){ $cpk_vol = -$cpk_vol }

		$vol_ratio = $cpk_vol/$dpk_vol;
		$atan_expr = sqrt(-$vol_ratio);
		$arctan = atan($atan_expr);
		$jcoup = (atan($atan_expr))/(2*$pi*$tau);
		if ($mult_tes eq "Y" or $mult_test eq "y"){
			$jcoup = $jcoup*1.1;
		}
		#print OUTFILE "expression: $atan_expr arctangent: $arctan J: $jcoup\n";
		print OUTFILE "$resname{$x} $x   $atm:HN   $jcoup\n";
	}
}

close OUTFILE;