frequency: Process multiple reflections, do plotting in Perl script rather than shell

1 files changed, 136 insertions, 20 deletions

diff --git a/scripts/frequency b/scripts/frequency
index 1f9c66e6..0b2feb56 100755
--- a/scripts/frequency
+++ b/scripts/frequency
@@ -3,49 +3,165 @@
 use strict;
 
 my $bins = 200;
+my $nplots = 32;
+
 open(FH, $ARGV[0]);
 
-my $min = 99999999999999999999.0;
-my $max = 0.0;
 my $line;
+my %counts;
+
+# --- Step 1: Count reflections and find the most popular ones
 
 while ( $line = <FH> ) {
 
+	my $h;
+	my $k;
+	my $l;
+
 	chomp($line);
-	if ( $line < $min ) {
-		$min = $line;
-	}
-	if ( $line > $max ) {
-		$max = $line;
+
+	if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)/ ) {
+
+		my $key;
+
+		$h = $1;
+		$k = $2;
+		$l = $3;
+
+		$key = $h.",".$k.",".$l;
+		$counts{$key}++;
+
 	}
 
 }
 
-printf("%f -> %f\n", $min, $max);
+my $ref;
+my $i;
+$i = 0;
+foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts )
+{
+	$i++;
+	if ( $i > $nplots ) {
+		delete($counts{$ref});
+	} else {
+		printf("%s %i\n", $ref, $counts{$ref});
+	}
+}
 
-seek(FH, 0, 0);
 
-my $binsize = ($max+1-$min)/$bins;
-my @hist;
-my $i;
+# --- Step 2: Work out histogram bins for each reflection
 
-for ( $i=0; $i<$bins; $i++ ) {
-	$hist[$i] = 0;
-}
+seek(FH, 0, 0);
+
+my %mins;
+my %maxs;
 
 while ( $line = <FH> ) {
 
-	my $bin;
+	my $key;
+	my $intensity;
 
 	chomp($line);
 
-	$bin = int(($line-$min)/$binsize);
-	$hist[$bin]++;
+	if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)\s+([\d\.]+)/ ) {
+
+		my $h;
+		my $k;
+		my $l;
+
+		$h = $1;
+		$k = $2;
+		$l = $3;
+		$intensity = $4;
+
+		$key = $h.",".$k.",".$l;
+
+		next unless exists($counts{$key});
+
+		unless ( exists($mins{$key}) ) {
+			$mins{$key} = $intensity;
+		}
+		unless ( exists($maxs{$key}) ) {
+			$maxs{$key} = $intensity;
+		}
+		if ( $intensity < $mins{$key} ) {
+			$mins{$key} = $intensity;
+		}
+		if ( $intensity > $maxs{$key} ) {
+			$maxs{$key} = $intensity;
+		}
+
+	}
 
 }
 
-for ( $i=0; $i<$bins; $i++ ) {
-	printf("%f\t%i\n", $min+(($i+0.5)*$binsize), $hist[$i]);
+# --- Step 3: Bin the counts and plot
+
+open(GP, "| gnuplot");
+printf(GP "set term postscript enhanced font \"Helvetica,20\"\n");
+printf(GP "set output \"histo.ps\"\n");
+foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts )
+{
+	my $max = $maxs{$ref};
+	my $min = $mins{$ref};
+	my $binsize = ($max+1-$min)/$bins;
+	my @hist;
+
+	printf("%s %f -> %f\n", $ref, $mins{$ref}, $maxs{$ref});
+
+	for ( $i=0; $i<$bins; $i++ ) {
+		$hist[$i] = 0;
+	}
+
+	seek(FH, 0, 0);
+	while ( $line = <FH> ) {
+
+		my $bin;
+		my $h;
+		my $k;
+		my $l;
+		my $intensity;
+
+		chomp($line);
+
+		if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)\s+([\d\.]+)/ ) {
+
+			my $key;
+
+			$h = $1;
+			$k = $2;
+			$l = $3;
+			$intensity = $4;
+
+			$key = $h.",".$k.",".$l;
+
+			next unless ( $ref eq $key );
+
+			$bin = int(($intensity-$min)/$binsize);
+			$hist[$bin]++;
+
+		}
+
+	}
+
+	open(OFH, "> ".$ref.".dat");
+	for ( $i=0; $i<$bins; $i++ ) {
+		printf(OFH "%f\t%i\n", $min+(($i+0.5)*$binsize), $hist[$i]);
+	}
+	close(OFH);
+
+	my $nref = $ref;
+	$nref =~ s/\,/\ /g;
+	my $title = "t \"".$nref." reflection\"";
+	printf(GP "plot \"".$ref.".dat\" u 1:2 w histeps ".$title."\n");
+
 }
 
 close(FH);
+close(GP);
+
+foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts ) {
+	unlink($ref.".dat");
+}
+system("ps2pdf histo.ps");
+unlink("histo.ps");