aboutsummaryrefslogtreecommitdiff
path: root/scripts/frequency
diff options
context:
space:
mode:
Diffstat (limited to 'scripts/frequency')
-rwxr-xr-xscripts/frequency156
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");