#!/usr/bin/perl

# USAGE makeIntegral peak.list 2d.txt 

$scale=1;
$field = 500;  #nmr field MHz for scaling Hz position change

# NOTHING TO EDIT BELOW HERE:
use Math::Trig;
if($#ARGV < 0)
{
   $prgName = `basename $0`;
   chop($prgName);
   die <<EOHELP;

Program reads in 2d.txt file from kinetics expts. ie 1Ds separated by an "&"
and creates intensity, max position and integrals from a peak.list file. 
files created: intensity.txt, integral.txt, position.txt, Int

peak.list format (count ppm_range):
0	1.3:2.4
1	4.5:5.2
...

2d.txt format (ppm intensity):
1.223	1002.11
1.224	1002.00
...

Usage: $prgName peak.list 2d.txt 
EOHELP
}

open LIST, $ARGV[0];
$i = 0;
while (defined($line = <LIST>)) {;
	@fields = split(/\s+/, $line);
			$peak[$i] = $fields[0];
			($a[$i], $b[$i]) = split(/:/, $fields[1]);
    			if ($a==$b){                     
			$high[$i]=$a[$i];
			$low[$i]=$b[$i];}
    			if ($a[$i]>$b[$i]){                     
			$high[$i]=$a[$i];
			$low[$i]=$b[$i];}
			else{
			$high[$i]=$b[$i];
			$low[$i]=$a[$i];
			}
                        $high[$i] = sprintf("%.3f", $high[$i]);
                        $low[$i] = sprintf("%.3f", $low[$i]);
#			print "i is $i, peak $peak[$i], ppm is $high[$i] $low[$i]\n";
			$i += 1;
			$numpks=$i;
};
#			print "number of pks $numpks\n";
close LIST;

#OPEN FILE 2; 2d.txt
$file = $ARGV[1];
$count=0;
$spectrum=1;
open(IN, $file) || die "Can't open file $file.\n";
while(<IN>)
{
   next if /^#/;
   @value = split;

   $ppm{$spectrum}{$count} = $value[0];
   $int{$spectrum}{$count} = $value[1];

   $numList{$spectrum}{$count} = $int;

   if($value[1] > $max{$spectrum}){$max{$spectrum} = $value[1];}
   if($count==0){$max{$spectrum} = $value[1];}
   $count++;
   if($value[0] eq '&'){ $spectrum++; $cnt_max=$count-1; $count=0; $spec_cnt++}
   #printf "%4.3f %4.3f  %4.3f \n ",$value[0], $value[1], $spectrum;

}
close(IN);

## sort the data, integrate and determine maximums.
for ($i=0; $i < $numpks; $i++){
foreach $spectrum ( sort numerically keys %ppm ) {
foreach $count ( sort numerically keys %{$ppm{$spectrum} } ) {
    if ($ppm{$spectrum}{$count} < $high[$i] && $ppm{$spectrum}{$count}> $low[$i] && $spectrum > 0){
	    #printf "%7.3f \t", $ppm{$spectrum}{$count};
	    #printf "%7.3f \n", $int{$spectrum}{$count};
       $add[$i][$spectrum] =$add[$i][$spectrum] +  $int{$spectrum}{$count};	#integrate the peak regions
       #if($int{$spectrum}{$count} > $maxint[$i][$spectrum]){$maxint[$i][$spectrum] = $int{$spectrum}{$count} ; $maxpos[$i][$spectrum] = $ppm{$spectrum}{$count};}
       if(abs($int{$spectrum}{$count}) > abs($maxint[$i][$spectrum])){$maxint[$i][$spectrum] = $int{$spectrum}{$count} ; $maxpos[$i][$spectrum] = $ppm{$spectrum}{$count};}
		}
	}
   }
}

######################################
# PRINT OUT NORMALIZED INTEGRATED VALUES to INT file for glove
$numfiles=@spectrum;
 open(OUTPUT, ">Int");
for ($i=0; $i < $numpks; $i++){
       printf OUTPUT " %-2.0f %-2.0f HN ",$i, $i;
for ($x=1; $x<= $spec_cnt; $x++){
	printf OUTPUT "%2.5f ", $add[$i][$x]/abs($add[$i][1]);	#print normalized to 1st
	#printf OUTPUT "%2.5f ", $add[$i][$x];			#print absolute values
	}
       print OUTPUT "\n";
}
close(OUTPUT);

######################################
# PRINT OUT INTEGRATED VALUES to integrals.txt file
$numfiles=@spectrum;
open(OUTPUT2, ">integral.txt");
printf OUTPUT2 "#X	Y0	Y1	etc\n";
for ($x=1; $x<= $spec_cnt; $x++){
       printf OUTPUT2 "%3.0f\t",$x;
for ($i=0; $i < $numpks; $i++){
	printf OUTPUT2 "%2.5f\t", $add[$i][$x]/abs($add[$i][1]);
	#printf OUTPUT2 "%2.5f\t", $add[$i][$x];
	}
       print OUTPUT2 "\n";
}
close(OUTPUT2);

## PRINT OUT intensity.txt file
$numfiles=@spectrum;
open(OUTPUT3, ">intensity.txt");
printf OUTPUT3 "#X	Y0	Y1	etc\n";
for ($x=1; $x<= $spec_cnt; $x++){
       printf OUTPUT3 "%3.0f\t",$x;
for ($i=0; $i < $numpks; $i++){
	printf OUTPUT3 "%2.5f\t", $maxint[$i][$x]/abs($maxint[$i][1]);
       #printf OUTPUT3 "%2.5f\t", $maxint[$i][$x];
	}
       print OUTPUT3 "\n";
}
close(OUTPUT3);

# PRINT OUT Max position position.txt file
$numfiles=@spectrum;
open(OUTPUT4, ">position.txt");
printf OUTPUT4 "#X	X0	X1	etc\n";
for ($x=1; $x<= $spec_cnt; $x++){
       printf OUTPUT4 "%3.0f\t",$x;
for ($i=0; $i < $numpks; $i++){
       #printf OUTPUT4 "%2.5f\t", $maxpos[$i][$x];
       #$pos= ($maxpos[$i][$x] - $maxpos[$i][1])*$field;
       printf OUTPUT4 "%2.5f\t", $maxpos[$i][$x];
       #printf OUTPUT4 "%2.5f\t", $pos;
	}
       print OUTPUT4 "\n";
}
close(OUTPUT4);

exit;
sub numerically { $a <=> $b; }