#!/usr/bin/perl
$numres = 21;
$done = 0;

while ($line = <STDIN>)  {
    chop($line);
    @volinfo = split(' ', $line); 
    $vv = @volinfo;
    $resfile = $ARGV[0]; 
    
    if ($resfile && !($done))  {
	&read_in_types;
	&create_Vstats;
	$done++;
    } elsif (!($resfile) && ($volinfo[0] eq "#"))  {
	$ll = $vv - 1;
	$resfile = $volinfo[$ll];
	&read_in_types;
	&create_Vstats;
    } elsif($volinfo[0] ne "" && $volinfo[0] ne "Total" && $volinfo[0] ne "atom")  {
	if($vv == 1)  {
	    $res = $volinfo[0];
	} else  {
	    $resdef = "$res $volinfo[0]";
	    $typ = $atomdefs{$resdef};
	    # for the sum mean and deviation
	    $Vnum{$typ} += $volinfo[1];
	    $Vtot{$typ} += $volinfo[1]*$volinfo[2];
	    $Vstd{$typ} += $volinfo[1]*($volinfo[3]**2+$volinfo[2]**2);
	    if($volinfo[5] < $Vmin{$typ})  { $Vmin{$typ} = $volinfo[5]; }
	    if($volinfo[6] > $Vmax{$typ})  { $Vmax{$typ} = $volinfo[6]; }
	    # for mean and deviation from outfile
	    $Vnn{$typ} += 1;
	    $Vsu{$typ} += $volinfo[2];
	    $Vsq{$typ} += $volinfo[2]**2;
	    if($volinfo[2] < $Vmi{$typ})  { $Vmi{$typ} = $volinfo[2]; }
	    if($volinfo[2] > $Vma{$typ})  { $Vma{$typ} = $volinfo[2]; }
	}
    }			
}

$sumVstd = 0.0; $cntVstd = 0.0;
$sumSstd = 0.0; $cntSstd = 0.0;
$count = 0;

print "#             ---------------- Sum ----------------    ---------- from outfile -----------\n";
print "# atomtype    num     mean    dev     min     max      num   mean     dev     min    max\n";
foreach $typ (sort keys(%Vnum))  {
    if(!($Vnum{$typ} == 0.0))  {
	$count++;
	# Sum from calculation -------------------------------------
	$vol = $Vtot{$typ}/$Vnum{$typ};
	$std = sqrt(($Vstd{$typ}/$Vnum{$typ})-($vol**2));
	$sumVstd += $Vnum{$typ}*($std**2); $cntVstd += $Vnum{$typ};
	printf "  %8s %8d %7.3lf %7.3lf %7.3lf %7.3lf", 
	$typ, $Vnum{$typ}, $vol, $std, $Vmin{$typ}, $Vmax{$typ}; 
	# from outfile ---------------------------------------------
	$vol = $Vsu{$typ}/$Vnn{$typ};
	$std = sqrt(($Vsq{$typ}/$Vnn{$typ})-($vol**2));
	$sumSstd += $Vnn{$typ}*($std**2); $cntSstd += $Vnn{$typ};
	printf "    %3d %7.3lf %7.3lf %7.3lf %7.3lf\n", 
	$Vnn{$typ}, $vol, $std, $Vmi{$typ}, $Vma{$typ}; 
    }
}

#getting resfile name alone 
@daname=split(/\//,$resfile);
$last = @daname;

$sVstd = sqrt($sumVstd/$cntVstd);
$sSstd = sqrt($sumSstd/$cntSstd);

printf "#        resfile  num        Sumdev                                 Sumdev\n";
printf "# Totals %8s %3d       %7.3lf                                %7.3lf\n", 
    $daname[$last-1], $count, $sVstd, $sSstd;


# residue file must be in a default format
# only doing the first 21 residues in file
sub read_in_types  {
    open(RES,$resfile) || die "couldn't open or find $resfile\n";
    print "# for residue typefile in $resfile\n";
    $jj=0;
    while($jj<$numres)  {
	$rline = <RES>;
	chop($rline);
	@res = split(' ', $rline);
	$rr = @res;
	if (  $rr > 0 )  {
	    for($ii=0;$ii<$res[1];$ii++)  {
		$rline = <RES>;
		@atoms = split(' ', $rline);
		if( !($atoms[1] =~ /=/) )  {
		    $resdef = "$res[0] $atoms[0]";
		    $atomdefs{$resdef} = $atoms[1];
		    $Vtyp{$atoms[1]} = 0.0;
# for old version of dealing with the cysteines
#		    if($res[0] eq "CYS" && $atoms[0] eq "SG")  {
#			$resdef = "$res[0] SGS";
#			$atomdefs{$resdef} = "S2H"; 
#			$Vtyp{"S2H"} = 0.0;
#		    }
		}
	    }
	    $jj++;
	}
    }				
    close(RES);	
}

# creates volume stats read in from stdin (from output of nVradii)
sub create_Vstats  {
    foreach $typ (sort keys(%Vtyp))  {
	# These find sum means and deviations of an atomtype from calc. over pdb
	$Vnum{$typ} = 0.0;
	$Vtot{$typ} = 0.0;
	$Vstd{$typ} = 0.0;
	$Vmin{$typ} = 1000000000;
	$Vmax{$typ} = -1000000000;
	# These find mean and deviation of atomtype from output file
	$Vsq{$typ} = 0.0;
	$Vnn{$typ} = 0.0;
	$Vsu{$typ} = 0.0;
	$Vmi{$typ} = 10000000000;
	$Vma{$typ} = -10000000000;
    }
}

