From 30012bfefa8c388c86b1fe0078fc3665798cfcc8 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Thu, 11 Nov 2010 17:11:45 +0100
Subject: Update scripts

---
 scripts/cell-please | 130 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 130 insertions(+)
 create mode 100755 scripts/cell-please

(limited to 'scripts/cell-please')

diff --git a/scripts/cell-please b/scripts/cell-please
new file mode 100755
index 00000000..bb8b70ce
--- /dev/null
+++ b/scripts/cell-please
@@ -0,0 +1,130 @@
+#!/usr/bin/perl -w
+
+use strict;
+use Fcntl;
+use POSIX;
+use IO::Handle;
+
+open(FH, $ARGV[0]);
+
+open(GP, "| gnuplot");
+autoflush GP 1;
+print(GP "set term postscript enhanced font \"Helvetica,20\"\n");
+print(GP "set output \"unitcell.ps\"\n");
+print(GP "unset key\n");
+print(GP "set xtics nomirror out rotate by 0\n");
+print(GP "unset xdata\n");
+print(GP "set format x \"% g\"\n");
+print(GP "set ylabel \"Frequency\"\n");
+print(GP "unset key\n");
+print(GP "set border lw 2\n");
+print(GP "set xtics nomirror out rotate by -60\n");
+print(GP "set grid\n");
+
+my $a;
+my $b;
+my $c;
+my $al;
+my $be;
+my $ga;
+
+print(GP "set xlabel \"Axis length / nm\"\n");
+print(GP "set title \"a\"\n");
+$a = &find_max(*FH, "^Cell\ parameters\ ([0-9\.]+)\ [0-9\.]+\ [0-9\.]+", "[0:60]");
+print(GP "set title \"b\"\n");
+$b = &find_max(*FH, "^Cell\ parameters\ [0-9\.]+\ ([0-9\.]+)\ [0-9\.]+", "[0:60]");
+print(GP "set title \"c\"\n");
+$c = &find_max(*FH, "^Cell\ parameters\ [0-9\.]+\ [0-9\.]+\ ([0-9\.]+)", "[0:60]");
+
+print(GP "set xlabel \"Angle / deg\"\n");
+print(GP "set title \"alpha\"\n");
+$al = &find_max(*FH, "([0-9\.]+)\ [0-9\.]+\ [0-9\.]+ deg\$", "[0:180]");
+print(GP "set title \"beta\"\n");
+$be = &find_max(*FH, "[0-9\.]+\ ([0-9\.]+)\ [0-9\.]+ deg\$", "[0:180]");
+print(GP "set title \"gamma\"\n");
+$ga = &find_max(*FH, "[0-9\.]+\ [0-9\.]+\ ([0-9\.]+) deg\$", "[0:180]");
+
+close(FH);
+close(GP);
+
+printf("Axis lengths: %5.2f %5.2f %5.2f nm\n", $a, $b, $c);
+printf("Angles: %5.2f %5.2f %5.2f deg\n", $al, $be, $ga);
+
+exit(0);
+
+
+
+sub find_max()
+{
+	my $fh = shift;
+	my $exp = shift;
+	my $lims = shift;
+
+	my $line;
+	my $hsteps = 1000;
+	my @hist;
+	my $hmin;
+	my $hmax;
+	my $first = 1;
+
+	seek $fh, 0, SEEK_SET or die("Couldn't rewind input");
+	while ( $line = <$fh> ) {
+
+		chomp $line;
+
+		if ( $line =~ /$exp/ ) {
+
+			my $val = $1;
+
+			if ( $first ) {
+				$hmin = $val;
+				$hmax = $val;
+				$first = 0;
+			}
+
+			if ( $val > $hmax ) { $hmax = $val; }
+			if ( $val < $hmin ) { $hmin = $val; }
+		}
+
+	}
+
+	my $hrange = $hmax - $hmin;
+
+	my $hstep = ($hmax - $hmin)/$hsteps;
+	for ( my $i=0; $i<$hsteps; $i++ ) {
+		$hist[$i] = 0;
+	}
+	seek $fh, 0, SEEK_SET or die("Couldn't rewind input");
+	while ( $line = <$fh> ) {
+		chomp $line;
+		if ( $line =~ /$exp/ ) {
+
+			my $val;
+			my $bin;
+
+			$val = $1;
+			$bin = floor(($val-$hmin)/$hstep);
+			$hist[$bin]++;
+
+		}
+	}
+	open(OFH, "> histogram.minated");
+	for ( my $i=0; $i<$hsteps; $i++ ) {
+		printf(OFH "%f %f\n", $hmin+$hstep*$i+($hstep/2), $hist[$i]);
+	}
+	close(OFH);
+	print(GP "plot ".$lims." [] \"histogram.minated\" u 1:2 w histeps lw 5 lc 3\n");
+
+	my $max = 0;
+	my $mval = 0;
+	for ( my $bin=0; $bin<$hsteps; $bin++ ) {
+		if ( $hist[$bin] > $mval ) {
+			$mval = $hist[$bin];
+			$max = $bin;
+		}
+	}
+	sleep(1);
+	unlink("histogram.minated");
+
+	$hmin + $hstep*$max;
+}
-- 
cgit v1.2.3