#!/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;
}