#!/usr/bin/perl -w 
use Getopt::Std;
use Text::Tabs;
require 'ctime.pl';
require "/home/vossman/lib/fixnum.pl";

getopt("ior", \%args);
$isresi{'CYT'}=1;
$isresi{'URI'}=1;
$isresi{'GUA'}=1;
$isresi{'ADE'}=1;
@times = split(' ', &ctime(time));
$time = $times[1].'-'.$times[2].'-'.$times[4];

#Set input variables
$in = $args{i} || "in.out";
$resi = $args{r} || "../../spec-defs.dat";   #designed for resi-defs.dat
if($args{o}) {
   $out = $args{o};
   open(OUT, "> $out")
      or die "Couldn't open $out for writing: $!\n";
} else {
   $out = "./stdvol258_$time.dat";
   open(OUT, "> $out")
      or die "Couldn't open $out for writing: $!\n";
}
print "Writing to $out...\n";

#$header = " atom\ttype\tcount\t\% OK\t mode\t mean\t sd\t% SD\t  min\t  max\n";
$header = "";
#$header2 = "------\t-----\t-----\t----\t-----\t------\t----\t----\t------\t------\n";
$header2 = "";
print OUT 
"# NucProt type set of 258 Volumes\n#\n# Generated: ",&ctime(time)
,"# Input file $in\n# Resi file $resi\n#\n\n";

print OUT "   atom  Count  Mean volume    SD       SD()    Minimum    Maximum\n";

#Open Files
open(RESI, "< $resi")
   or die "Couldn't open $resi for reading: $!\n";

$i=0;
while(defined($line = <RESI>)) {
   if($line && substr($line,0,1) && !$i) {
      $line =~ s/\s+$//; $line =~ s/^\s+//;
      @chunks = split(' ', $line);
      foreach $test (@chunks) { $test =~ s/\s+$//; $test =~ s/^\s+//; }
      if($chunks[1]) {
         push(@resi,$chunks[0]);
         $curresi = $chunks[0];
         $i=$chunks[1];
      }
   } elsif($curresi && $line && substr($line,2,1)) {
      $line =~ s/\s+$//; $line =~ s/^\s+//;
      @chunks = split(' ', $line);
      foreach $test (@chunks) { $test =~ s/\s+$//; $test =~ s/^\s+//; }
      push(@{$curresi},$chunks[0]);
      $i--;
   }
} close(RESI);

foreach $RESI (@resi) {
   $tvol = 0;
   $tvolsq = 0;
   $tokay = 0; $tposs = 0; $tbad = 0;
   if($out ne "stdout") { print STDERR "\n$RESI\t[ "; }
   $k = -1;
   print OUT "$RESI\n",$header,$header2;
   foreach $ATOM (@{$RESI}) {
      $k++;
      if($out ne 'stdout') { print STDERR "$ATOM "; }
      $okay = 0; $poss=0; $bad=0;
      $ovol = 0;
      $ovolsq = 0;
      $maxvol = 0;
      $minvol = 1500;
      open(IN, "< $in")
         or die "Couldn't open $in for reading: $!\n";
      while( defined($line = <IN>) ) {
        @chunks = split('\t',$line);
        foreach $test (@chunks) { $test =~ s/\s+$//; $test =~ s/^\s+//; }
        if($ATOM eq $chunks[4] && ( ($RESI eq "SUG" && $isresi{$chunks[3]})
            || $RESI eq $chunks[3])) {
           if($chunks[13] eq "\[OK\]") {
              $okay++;
              $vol = $chunks[10];
              $ovol += $vol;
              $ovolsq = $ovolsq + $vol*$vol;
              if($vol > $maxvol) { $maxvol = $vol; }
              if($vol < $minvol) { $minvol = $vol; }
           } 
           elsif($chunks[13] eq "possible") { $poss++; }
           elsif($chunks[13] eq "bad") { $bad++; }
        } 
     }
     close(IN);
     $tokay += $okay;
     $tbad += $bad;
     $tposs += $poss;
     if($okay > 1) {
        $ovolsq = sqrt( ($ovolsq - $ovol*$ovol/$okay) / ($okay - 1) );
        $ovol = $ovol / $okay;
        $tvol += $ovol;
        $tvolsq += $ovolsq*$ovolsq;
     } elsif ($okay == 1) {
        $ovolsq = 0;
     }
     if($okay > 1) {
#FIND MODE
        $line = "   $ATOM\t ".fixnum50($okay)."    ".
                     fixnum23($ovol)."     ".fixnum13($ovolsq)."    ".
		     fixnum23($ovolsq/$ovol*100)."     ".
                     fixnum23($minvol)."     ".fixnum23($maxvol)."\n";
        expand($line);
     
        print OUT $line;
     } else {
        print STDERR "\be ";
     }
  }
  if($out ne 'stdout') { print STDERR "]"; }
  if($tokay > 1) {
     $tvolsq = sqrt($tvolsq);
     print OUT "\n Total   ",fixnum50($tokay),"   ",
                fixnum33($tvol),"    ",fixnum23($tvolsq),"    ",
                fixnum23($tvolsq/$tvol*100),"\n\n";
  }
}
print STDERR "\n\nIf an \"e\" follows the character; no \[OK\] atoms were found.\n\n";