diff options
Diffstat (limited to 'scripts')
| -rwxr-xr-x | scripts/frequency | 156 | ||||
| -rwxr-xr-x | scripts/histogram | 15 |
2 files changed, 136 insertions, 35 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"); diff --git a/scripts/histogram b/scripts/histogram deleted file mode 100755 index 329b775e..00000000 --- a/scripts/histogram +++ /dev/null @@ -1,15 +0,0 @@ -#!/bin/sh - -./frequency.pl 103.hist > hist.dat - -gnuplot << EOF -set term postscript enhanced font "Helvetica,20" -set output "histo.ps" - -plot [-2000:60000] "hist.dat" u 1:2 w l t "103 reflection" -EOF - -ps2pdf histo.ps -rm -f histo.ps hist.dat - -evince histo.pdf |