#!/usr/bin/perl

## Copyright 2008 by James Rath and the University of Texas at Austin
## This code is available under a Creative Commons Attribution-ShareAlike license.  You
## are free to copy and share this as you see fit; however, you must attribute any copy
## to me.  You are also free to create derivative works; however, you must acknowledge
## the use of this original, and you must share your work under the same terms you
## received this copy.  For details, see http://creativecommons.org/licenses/by-sa/3.0/
## For permission for any other use, please contact me at ratjamm@alum.mit.edu

use POSIX qw(ceil);


## one of these days i'll write usage instructions...
@usage = ("blah\n",
		  "blah\n",
		  "blah\n");


## name on the command line tells us what file(s) to look in
# it'll also give us a title name for the EPS
if ($#ARGV+1) {
	($ARGV[0] eq "--help")  &&  die @usage;
	$titlePrefix = shift @ARGV;
}


# loop over the list of eigenvector files
# produce an EPS file for each one
open(VECLIST,"ls -1 ./$titlePrefix.?.node |");
while (<VECLIST>) {
	chop;
	s/.node$//;
	$title = $_;


## slurp up the nodes
open(NODES,"< $title.node");
# First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)>
# skip possible comment lines
$headerline = <NODES>;
while ($headerline =~ /^\s*#/ || $headerline =~ /^\s*$/) {
	$headerline = <NODES>;
	die "all comments in nodes file?!" if (!$headerline);
}
# read and check the information
($nnodes,$dim,$nattr,$bdymarkp) = ($headerline =~ /\s*(\d+)\s*(\d+)\s*(\d+)\s*(\d+)\s*/);
die "naughty nodes file\n" if (!$nnodes || $dim!=2 || $nattr!=1 || $bdymarkp!=1);
# preallocate space for data
$#indices = $nnodes - 1;
$#xs = $nnodes - 1;
$#ys = $nnodes - 1;
$#nodeData = $nnodes - 1;
$#bdymarks = $nnodes - 1;
# Remaining lines: <vertex #> <x> <y> [attributes] [boundary marker]
$nLines = 0;
while (<NODES>) {
	next if ((/^\s*#/) || (/^\s*$/));
	($indices[$nLines],$xs[$nLines],$ys[$nLines],$nodeData[$nLines],$bdymarks[$nLines]) = split(' ', $_);
	if ($bdymarks[$nLines]==1) { $nodeData[$nLines] = 0.0 }
	$nLines++;
}
die "didn't read all the nodes?!\n" if ($nLines!=$nnodes);
close(NODES);

# sort the data just in case the nodes weren't written in order
@revIndices = sort { $indicies[$a] <=> $indicies[$b] } 0..$#indices;
@xs = @xs[@revIndices];
@ys = @ys[@revIndices];
@nodeData = @nodeData[@revIndices];
@bdymarks = @bdymarks[@revIndices];

# get the range of the grid points to use in picking an aspect ratio
($xMin,$xMax) = &minmax(@xs);
($yMin,$yMax) = &minmax(@ys);
$xScale = $xMax - $xMin;
$yScale = $yMax - $yMin;
if ($xScale < $yScale) {
    $xAspect = $xScale/$yScale;
    $yAspect = 1.0;
} else { # $yScale < $xScale
    $xAspect = 1.0;
    $yAspect = $yScale/$xScale;
}

# get the range of the data to use in picking a color palette
($dataMin,$dataMax) = &minmax(@nodeData);
## hack: don't shift zero, and scale by the infinity norm
if ( -$dataMin<$dataMax ) {
	$dataMin = -$dataMax;
} else {
	$dataMax = -$dataMin;
}
## end hack
$dataScale = $dataMax - $dataMin;


## write the EPS header
open(EPS,"> $title.eps");
# EPS header
print EPS "%!PS-Adobe-3.0 EPSF-3.0\n";
# always make the boundingBox six inches by six inches
printf EPS "%%%%BoundingBox: 0 0 %d %d\n", ceil($xAspect * 432.0), ceil($yAspect * 432.0);
printf EPS "%%%%HiResBoundingBox: 0 0 %.5g %.5g\n", $xAspect * 432.0, $yAspect * 432.0;
print EPS "%%Creator: Jim Rath's node to EPS\n";
print EPS "%%Title: $title\n";
($sec,$min,$hour,$day,$month,$year,$dotw,$doty,$isdst) = localtime(time);
printf EPS "%%%%CreationDate: %.4d/%.2d/%.2d %.2d:%.2d:%.2d\n", $year+1900,$month+1,$day,$hour,$min,$sec;
print EPS "%%DocumentData: Clean7Bit\n";
print EPS "%%LanguageLevel: 1\n";
print EPS "%%Pages: 1\n";
print EPS "%%EndComments\n\n";

# Page header
print EPS "%%Page: 1 1\n";
print EPS "%%BeginPageSetup\n";
printf EPS "%.5g %.5g scale\n", $xAspect * 432.0/$xScale, $yAspect * 432.0/$yScale;
printf EPS "%.5g %.5g translate\n", -$xMin, -$yMin;
# h s b sethsbcolor
print EPS "/s {1.0 1.0 sethsbcolor} def\n";
print EPS "/n {newpath} def\n";
#print EPS "/m {moveto} def\n";
#print EPS "/l {lineto} def\n";
#print EPS "/r {rlineto} def\n";
#print EPS "/cf {closepath fill} def\n";
# x y w h rectfill;
#print EPS "/rf {rectfill} def\n";
print EPS "/t {moveto lineto lineto closepath fill} def\n";
print EPS "%%EndPageSetup\n\n";


## slurp up the elements and write 'em out
open(ELTS,"< $titlePrefix.ele");
# for checking, read the number of elements from the first line
# First line: <# of triangles> <nodes per triangle> <# of attributes>
$headerline = <ELTS>;
while ($headerline =~ /^\s*#/ || $headerline =~ /^\s*$/) {
	$headerline = <ELTS>;
	die "all comments in elements file?!" if (!$headerline);
}
($ntri,$nodepertri,$nattr) = ($headerline =~ /\s*(\d+)\s*(\d+)\s*(\d+)\s*/);
die "naughty elements file\n" if (!$ntri || $nodepertri!=3 || $nattr!=0);
# read the elements, processing one-by-one
$nLines = 0;
while (<ELTS>) {
	next if ((/^\s*#/) || (/^\s*$/));
# Remaining lines: <triangle #> <node> <node> <node> ... [attributes]
	($tri,$v1,$v2,$v3) = split(' ', $_);
# put a triangle the size of each element filled with the appropriate color
## FIXME: we should check whether we have a boundary node and use 0.0 for that nodeData
##        currently we assume it was written correctly above in the .node file
## FIXME: using level 3 postscript, we could define a shading pattern that's a linear
##        gradient color fill according to the node values; however, by the time we
##        hit a million elements, each element is only about one pixel in a 1000x1000 image
# h s b sethsbcolor  --  NB: 0==red==1 so we only go 5/6ths of the way around in hue
	$v1 = $v1-1;
	$v2 = $v2-1;
	$v3 = $v3-1;
	$hue = (5.0/18.0)*($nodeData[$v1]+$nodeData[$v2]+$nodeData[$v3]-3.0*$dataMin)/$dataScale;
#	if ($hue > 5.0/6.0) { $hue = 5.0/6.0 }
#	if ($hue < 0.0) { $hue = 0.0 }
	printf EPS "%.5g s\n", $hue;
	printf EPS "n %.5g %.5g %.5g %.5g %.5g %.5g t\n", $xs[$v3],$ys[$v3], $xs[$v2],$ys[$v2], $xs[$v1],$ys[$v1];
	$nLines++;
}
# check that we got 'em all
die "didn't read all the triangles?!\n" if ($nLines!=$ntri);
close(ELTS);

## write the EPS footer and close up shop
# NB: good for broken interpreters when displaying pure EPS
print EPS "\nshowpage\n\n";
# Page/EPS footer
print EPS "%%PageTrailer\n";
print EPS "%%Trailer\n";
print EPS "%%EOF\n";
close(EPS);


}
# end while(<VECLIST>)
close(VECLIST);


### helper functions

sub minmax {
	local($min) = pop(@_);
	local($max) = $min;
	foreach $foo (@_) {
		$min = $foo if $min > $foo;
		$max = $foo if $max < $foo;
	}
	($min,$max);
}