#!/usr/bin/perl -nl

# USAGE: call it with all C files on the command line
# $ ./terra_export.pl c<casenum>.*.<timestep>
#
# Author : Christoph Moder <moder@geophysik.uni-muenchen.de>
# Version: 2011-10-19

use constant PI => atan2(0, -1);
use constant ROUND => 1;	# project points of subdivision on sphere or leave them on the icosahedron
use POSIX qw( _exit );

sub acos { return atan2( sqrt(1 - $_[0] * $_[0]), $_[0] ); }

# converts radians to degrees
sub rad2deg { $_[0] * 180 / atan2(0, -1) }

# limits a number to the range between min and max
sub limit {
	my ($val, $min, $max) = @_;

	($val < $min ? $min : $val > $max ? $max : $val );
}

# read input, remove comments (starting with "#" character), repeat until non-empty entry (or default value)
sub get_input {
	my ($prompt, $default, $allowed) = @_;

	READINPUT:
	{
		print STDERR "$prompt [Enter = $default]";

		$_ = <STDIN>;

		# replace empty input by default; remove comment and space
		/^$/ ? $_ = $default : (s/^\s*//, s/#.*//, s/\s*$//);
	 
		# go back if empty input (e.g. comment line in script)
		$_ eq "" and redo READINPUT;

		# if input in allowed range:
		# * if not interactively, print it to generate script
		# * otherwise do nothing
		# * return value
		# if not allowed:
		# * if interactively, then retry
		# * otherwise: output error message, switch to interactive mode and retry
		$_ =~ $allowed ? (-t STDOUT ? 0 : print( "$_\t# $ prompt"), $_) : -t STDIN ? redo READINPUT : ( print( STDERR "SCRIPT ERROR: terminating script, switching to manual input" ), close( STDIN ), open( STDIN, "/dev/tty" ) );
	}
}

# converts depth (in km or percent) to layer number; requires @layers
sub d2l {
	my $input = shift;
	my $i = 0;
	my $depth = ( $input =~ /%/ ? $input / 100 * 6370E3 : $input =~ /km/ ? 6370E3 - $_ * 1E3 : 0 );

	# if given in km or percent of radius:
	# 1. assigns indices to the layers (=> pairs: [index, depth])
	# 2. sorts them according to depth difference to desired depth
	# 3. removes depth, so only layer numbers are left
	# 4. takes first element => the layer number where the depth difference is minimal
	my $layer = ( $input =~ /%|km/ ? (map { $_->[0] } sort { $a->[1] <=> $b->[1] } map { [ $i++, abs($depth - $_)] } @layers)[0] : (0 + $input));

	return $layer;
}

# generates TERRA grid
sub grdgen {
	my $a   = 2 * acos( 1 / (2 * sin( PI / 5 ) ) );
	my $lvt = int( 1.45 * log( $nt ) );
        my @ip = ( -1, 1, 3, 5, -3, 0, 2, 4, -4, -2 );

	my $nd = 10;

	for my $diamond (0..$nd - 1)
	{
		my $sgn = (4.5 - $diamond <=> 0);	# upper or lower hemisphere
		my $phi = $ip[$diamond] * PI / 5;	# longitude shift

		# coordinates of the four corners of the diamond
		# first corner: pole (e.g. N)
		# second corner: south east of N pole
		# third corner: south west of N pole
		# fourth corner: south of N pole
		$xg[$diamond][0][0][0] = 0;
		$xg[$diamond][0][0][1] = 0;
		$xg[$diamond][0][0][2] = $sgn;

		$xg[$diamond][$nt][0][0] = sin( $a ) * cos( $phi + 2 * PI / 5 );
		$xg[$diamond][$nt][0][1] = sin( $a ) * sin( $phi + 2 * PI / 5 );
		$xg[$diamond][$nt][0][2] = cos( $a ) * $sgn;

		$xg[$diamond][0][$nt][0] = sin( $a ) * cos( $phi );
		$xg[$diamond][0][$nt][1] = sin( $a ) * sin( $phi );
		$xg[$diamond][0][$nt][2] = cos( $a ) * $sgn;

		$xg[$diamond][$nt][$nt][0] = sin( $a ) * cos( $phi + PI / 5 );
		$xg[$diamond][$nt][$nt][1] = sin( $a ) * sin( $phi + PI / 5 );
		$xg[$diamond][$nt][$nt][2] = -cos( $a ) * $sgn;

		# the diamond is iteratively subdivided
		# i1 direction: from pole to point xv[1] = row
		# i2 direction: from pole to point xv[2] = column
		# diagonal: between xv[1] and xv[2] = horizontal diagonal
		for my $k (0..$lvt - 1)
		{
			my $m = 2**$k;
			my $l = $nt / $m;

			# rows of diamond
			for my $j1 (0..$m)
			{
				for my $j2 (0..$m - 1)
				{
					my $i1 = $j1 * $l;
					my $i2 = $j2 * $l + $l / 2;

					# for x, y, z coordinates: interpolate from previous grid resolution
					$xg[$diamond][$i1][$i2][$_] = ($xg[$diamond][$i1][$i2 - $l / 2][$_] + $xg[$diamond][$i1][$i2 + $l / 2][$_]) / 2 for( 0..2 );

					if( ROUND )
					{
						# compute radius of new grid point
						my $xnorm = 1 / sqrt( $xg[$diamond][$i1][$i2][0]**2 + $xg[$diamond][$i1][$i2][1]**2 + $xg[$diamond][$i1][$i2][2]**2 );

						# project back onto sphere
						$xg[$diamond][$i1][$i2][$_] *= $xnorm for( 0..2 );
					}
				}
			}
			
			# columns of diamond
			for my $j1 (0..$m)
			{
				for my $j2 (0..$m - 1)
				{
					my $i1 = $j2 * $l + $l / 2;
					my $i2 = $j1 * $l;

					# for x, y, z coordinates: interpolate from previous grid resolution
					$xg[$diamond][$i1][$i2][$_] = ($xg[$diamond][$i1 - $l / 2][$i2][$_] + $xg[$diamond][$i1 + $l / 2][$i2][$_]) / 2 for( 0..2 );

					if( ROUND )
					{
						# compute radius of new grid point
						my $xnorm = 1 / sqrt( $xg[$diamond][$i1][$i2][0]**2 + $xg[$diamond][$i1][$i2][1]**2 + $xg[$diamond][$i1][$i2][2]**2 );
						
						# project back onto sphere
						$xg[$diamond][$i1][$i2][$_] *= $xnorm for( 0..2 );
					}
				}
			}
			
			# diagonals of diamond
			for my $j1 (0..$m - 1)
			{
				for my $j2 (0..$m - 1)
				{
					my $i1 = $j1 * $l + $l / 2;
					my $i2 = $j2 * $l + $l / 2;
					
					# for x, y, z coordinates: interpolate from previous grid resolution
					$xg[$diamond][$i1][$i2][$_] = ($xg[$diamond][$i1 - $l / 2][$i2 + $l / 2][$_] + $xg[$diamond][$i1 + $l / 2][$i2 - $l / 2][$_]) / 2 for( 0..2 );
					
					if( ROUND )
					{
						# compute radius of new grid point
						my $xnorm = 1 / sqrt( $xg[$diamond][$i1][$i2][0]**2 + $xg[$diamond][$i1][$i2][1]**2 + $xg[$diamond][$i1][$i2][2]**2 );
						
						# project back onto sphere
						$xg[$diamond][$i1][$i2][$_] *= $xnorm for( 0..2 );
					}
				}
			}
		}
	}
}

sub get_index
{
	my ($dataset, $layer, $diamond, $diamond_row, $diamond_col, $component) = @_;

	return( $grd_points_per_dataset * $dataset_component_offset[$dataset] + $layer * $grd_points_per_layer * $datavectorsize[$dataset] + $diamond * $grd_points_per_diamond * $datavectorsize[$dataset] + $diamond_row * $grd_points_per_row * $datavectorsize[$dataset] + $diamond_col * $datavectorsize[$dataset] + $component );
}


###########################################################################
# MAIN PROGRAM
###########################################################################

# not BEGIN, since the DATA handle is not available yet
INIT {
	# no parameters given, or "-h" => show documentation
	exec "perldoc $0" if( 0 == @ARGV or ( !-s $ARGV[0] and $ARGV[0] =~ /^-?-h/ ) );
	
	# one parameter "-e": print example file
	if( !-s $ARGV[0] and $ARGV[0] =~ /^-e$/ )
	{
		local $\;
		$|++;	# buffering must be switched off when using _exit, otherwise the output is lost

		print STDOUT while( ($_ = <DATA>) !~ /^__END__$/ );
		_exit(0);	# without END block
	}

	$, = "\t";
	@datavectorsize = (1, 3, 1, 1);	# temperature, velocity, pressure
	@datavectorname = ("temperature", "velocity", "pressure", "temp_deviation");

	$max_dataset = limit( get_input( "How many datasets should be read (@datavectorname[0..$#datavectorname-1])?", 3, qr/[1-3]/ ), 1, 3 );
	print STDERR "Your entry: " . ($max_dataset) . " => datasets: (@datavectorname[0..$max_dataset-1])";
	print STDERR "";

	print STDERR "NT", "NR", "CPUs", "SUBDIV", "D/CPU", "part.", "s/layer", "all/layer", "s", "all";
	print STDERR "=========================================================================================";

	# initializations
	$layer_lines = 100;
	$f = 0;		# file number (0 = we are in the first input file)
	$dataset = 0;
	my $i;
	@dataset_component_offset = map { ($i += $_) - $_ } @datavectorsize;
}


###########################################################################
# READING INPUT DATA
###########################################################################

# if input is compressed: close file and open via pipe from zcat
($. == 1) and /^\x1f\x8b/ and (close( ARGV ), open( ARGV, "zcat '$ARGV' |" ), next);

# first line of first file: read grid resolution; together with number of files on commandline (= number of processes) compute total resolution
if($f == 0 and $. == 1) {
	@F = split;

	$nr = $F[0] + 1;	# number of layers: from file header
	$cpus = 1 + @ARGV;	# number of processes: from number of files
	$diamonds_per_cpu = 10 - 5 * ((log($cpus)/log(2)) & 1);
	$subdivisions = sqrt( $cpus * $diamonds_per_cpu / 10 );	# number of partitions per diamond side
	$nt = $F[1] * $subdivisions;	# number of points per diamond side = length of partition * number of partitions per diamond side
	$layer_lines = int( $nr / 10 + 1 );	# length of layer depth entries in the header
	$header_lines = (1 + 4 + $layer_lines + 2);	# total length of file header

	# these dimensions are for the input files; in the final grid, the doubled nodes at the partition boundaries are omitted (except at diamond boundaries)
	$points_per_partition_side = $nt / $subdivisions + 1;
	$points_per_partition = $points_per_partition_side**2;
	$points_per_layer_per_cpu = $points_per_partition * $diamonds_per_cpu;
	$points_per_layer = $points_per_layer_per_cpu * $cpus;
	$points_per_cpu = $points_per_layer_per_cpu * $nr;
	$points = $points_per_layer * $nr;
	$dataset_line_offsets[$_] = ($_ ? $dataset_line_offsets[$_ - 1] + $header_lines + int( $points_per_cpu * $datavectorsize[$_ - 1] / 15 + 0.5 ) : 0) for( 0..@datavectorsize );

	$grd_points_per_row = $nt + 1;
	$grd_points_per_diamond = $grd_points_per_row**2;
	$grd_points_per_layer = $grd_points_per_diamond * 10;
	$grd_points_per_dataset = $grd_points_per_layer * ($nr + 1);
	
	# reserve memory at once by storing the last array element
	$data[ get_index( $max_dataset, $nr + 1, 9, $nt + 1, $nt + 1, $datavectorsize[$max_dataset] - 1 ) ] = 0;

	print STDERR $nt, $nr, $cpus, "${subdivisions}x${subdivisions}", $diamonds_per_cpu, $points_per_partition, $points_per_layer_per_cpu, $points_per_layer, $points_per_cpu, $points;

	print STDERR "generating grid ...";
	grdgen();

	print STDERR "reading data files ...";
}

# store depth layers (only in first input file)
($f == 0 and $. > 5 and $. <= (5 + $layer_lines) ) and @F = split and @layers = (@layers, @F);

# read data
if( $. > $header_lines + $dataset_line_offsets[$dataset] and $. <= $dataset_line_offsets[$dataset + 1] )
{
	# file format: 15 fields, each 10 characters wide
	@F = unpack( "(A10)15", $_ );

	for( @F )
	{
		$partition_in_diamond = $f;	# rectangular partition number (= MPI process), same on different diamonds
		$diamond = int( $fileoffset / $points_per_partition ) % $diamonds_per_cpu + 5 * int( $partition_in_diamond / $subdivisions**2 );

		$partition_in_diamond_row = int( $partition_in_diamond / $subdivisions ) % $subdivisions;
		$partition_in_diamond_col = $partition_in_diamond % $subdivisions;

		$partition_offset = int( $fileoffset ) % $points_per_partition;	# offset within the current partition (e.g. 1 = second number in the current partition)

		$partition_row = int( $partition_offset / $points_per_partition_side );	# row within the current partition
		$partition_col = $partition_offset % $points_per_partition_side;	# column within the current partition

		$layer = int( $fileoffset / ($diamonds_per_cpu * $points_per_partition * $datavectorsize[$dataset]) ) % (0 + @layers);	# radial layer; 0 = surface

		$component = int( $fileoffset / ($diamonds_per_cpu * $points_per_partition) ) % $datavectorsize[$dataset];

		$diamond_row = $partition_in_diamond_row * ($points_per_partition_side - 1) + $partition_row;
		$diamond_col = $partition_in_diamond_col * ($points_per_partition_side - 1) + $partition_col;

		# skip doubled points at partition boundaries; but not at diamond boundaries
		$do_skip = ( ($partition_row == $points_per_partition_side - 1 and $partition_in_diamond_row < $subdivisions - 1) or ($partition_col == $points_per_partition_side - 1 and $partition_in_diamond_col < $subdivisions - 1) ? 1 : 0 );

		# save data values	
		$data[ get_index( $dataset, $layer, $diamond, $diamond_row, $diamond_col, $component ) ] = $_ unless $do_skip;

		# record mean temperature
		$mean[$layer] += $_ / $points_per_layer unless( $do_skip or $dataset != 0 );

		$fileoffset++;

		# switch to next dataset if all data points are read
		($dataset++, $fileoffset = 0 ) if( $fileoffset >= $points_per_cpu * $datavectorsize[$dataset] );
		close( ARGV ) if( $dataset > $max_dataset );
	}
}

# reset line counter and increase file counter ($f) at end of a file
eof and (($. = ++$f - $f), $fileoffset = 0, $dataset = 0, print STDERR "reading file $f");


END {
###########################################################################
# USER INTERFACE LOOP
###########################################################################
	
	# default values
	$output_filenumber = 0;
	$input_fileformat  = "p";
	$input_diamonds    = "0-9";
	$input_horiz_x     = "0-$nt";
	$input_horiz_y     = "0-$nt";
	$input_layers      = "0-" . ($nr - 1);
	$input_downsample  = 1;
	$input_datasets    = "1,0,0,1";
	$input_filename    = "outputfile%i_%d_%l";

	while(1)
	{
		my @input;

		# FILE FORMAT
		$input_fileformat = get_input( "File format: (p)araview, (g)ocad, (x)yz? Or (q)uit program?", $input_fileformat, qr/^[pgxq]$/i );
		exit if( $input_fileformat =~ /q/i );


		# DIAMONDS
		$input_diamonds = get_input( "Which diamonds do you want to export (0-9; e.g. \"2-5\" or \"1,3,8-9\")?", $input_diamonds, qr/[0-9, -]+/ );
		@input = split /,/, $input_diamonds;
		my @write_diamonds = ();
		for( @input )
		{
			my @b;
			@b = split( /(-)/ ), (@b >= 2) and ( push( @write_diamonds, ( limit( 0 + $b[0], 0, 9 ) .. limit( ($b[2] ne '' ? 0 + $b[2] : 9 ), 0, 9 ) ) ), next );
			push( @write_diamonds, ( limit( 0 + $_, 0, 9 ) ) );
		}
		# remove duplicate entries
		@write_diamonds = sort { $a <=> $b } keys %{ { map { $_ => 1 } @write_diamonds } };
		print STDERR "Your entry: @write_diamonds";


                # HORIZONTAL RESOLUTION
		$input_horiz_x = get_input( "Which lines in x direction?", $input_horiz_x, qr/[0-9,+ -]+/ );
		@input = split /,/, $input_horiz_x;
		my @write_horiz_x = ();
		for( @input )
		{
			my @b;
			@b = split( /(-)/ ), (@b >= 2) and ( push( @write_horiz_x, ( limit( 0 + $b[0], 0, $nt ) .. limit( ($b[2] ne '' ? 0 + $b[2] : $nt ), 0, $nt ) ) ), next );
			push( @write_horiz_x, ( limit( 0 + $_, 0, $nt ) ) );
		}
		# remove duplicate entries
		@write_horiz_x = sort { $a <=> $b } keys %{ { map { $_ => 1 } @write_horiz_x } };

		$input_horiz_y = get_input( "Which lines in y direction?", $input_horiz_y, qr/[0-9,+ -]+/ );
		@input = split /,/, $input_horiz_y;
		my @write_horiz_y = ();
		for( @input )
		{
			my @b;
			@b = split( /(-)/ ), (@b >= 2) and ( push( @write_horiz_y, ( limit( 0 + $b[0], 0, $nt ) .. limit( ($b[2] ne '' ? 0 + $b[2] : $nt ), 0, $nt ) ) ), next );
			push( @write_horiz_y, ( limit( 0 + $_, 0, $nt ) ) );
		}
		# remove duplicate entries
		@write_horiz_y = sort { $a <=> $b } keys %{ { map { $_ => 1 } @write_horiz_y } };

		print STDERR "Your entry (x): @write_horiz_x";
		print STDERR "Your entry (y): @write_horiz_y";


		# LAYERS
		$input_layers = get_input( "Which layers do you want to export (e.g. \"1-7\", \"100km + 3\", \"95% - 90%\")?", $input_layers, qr/[0-9,+%km -]+/ );
		@input = split /,/, $input_layers;
		my @write_layers = ();
			
		# process comma separated list
		for( @input )
		{
			my @b;

			# if there are plus signs in the current element: push range from first element to first plus second element
			@b = split( /(\+)/ ), (@b >= 2) and ( push( @write_layers, ( limit( d2l( $b[0] ), 0, @layers - 1 ) .. limit( d2l($b[0]) + ($b[2] ne '' ? $b[2] : 1), 0, @layers ) ) ), next );

			# if there are minus signs: similar, but consider second part as absolute
			@b = split( /(-)/ ), (@b >= 2) and ( push( @write_layers, ( limit( d2l( ($b[0] ne '' ? $b[0] : 0) ), 0, @layers - 1 ) .. limit( d2l( ($b[2] ne '' ? $b[2] : @layers - 1) ), 0, @layers - 1 ) ) ), next );
			push( @write_layers, ( limit( d2l($_), 0, @layers - 1 ) ) );
		}

		# remove duplicate entries
		@write_layers = sort { $a <=> $b } keys %{ { map { $_ => 1 } @write_layers } };

		# only one layer: append or prepend the next layer
		#push( @write_layers, $write_layers[-1] + 1 ) if( 1 == @write_layers and $write_layers[-1] < @layers - 1 );
		#unshift( @write_layers, $write_layers[-1] - 1 ) if( 1 == @write_layers and $write_layers[-1] == @layers - 1 );

		print STDERR "Your entry: @write_layers";
		

		# DOWNSAMPLING
		$input_downsample = limit( get_input( "Downsampling factor?", $input_downsample, qr/[1-9][0-9]*/ ), 1, @{[sort { $a <=> $b } ( 0 + @write_layers, 0 + @write_horiz_x, 0 + @write_horiz_y )]}[-1] );
		print STDERR "your input: $input_downsample; sort: ", @{[sort { $a <=> $b } ( 0 + @write_layers, 0 + @write_horiz_x, 0 + @write_horiz_y )]}[0];
		# downsample
		@write_layers = grep { ($_ - $write_layers[0]) % $input_downsample == 0 } @write_layers;
		@write_horiz_x = grep { ($_ - $write_horiz_x[0]) % $input_downsample == 0 } @write_horiz_x;
		@write_horiz_y = grep { ($_ - $write_horiz_y[0]) % $input_downsample == 0 } @write_horiz_y;
		# only one layer: append or prepend the next layer
		( 1 == @write_layers ? $write_layers[-1] < @layers - 1 ? push( @write_layers, $write_layers[-1] + 1 ) : unshift( @write_layers, $write_layers[-1] - 1 ) : 0 );
		( 1 == @write_horiz_x ? $write_horiz_x[-1] < @layers - 1 ? push( @write_horiz_x, $write_horiz_x[-1] + 1 ) : unshift( @write_horiz_x, $write_horiz_x[-1] - 1 ) : 0 );
		( 1 == @write_horiz_y ? $write_horiz_y[-1] < @layers - 1 ? push( @write_horiz_y, $write_horiz_y[-1] + 1 ) : unshift( @write_horiz_y, $write_horiz_y[-1] - 1 ) : 0 );
		print STDERR "Your entry (x): @write_horiz_x\nYour entry (y): @write_horiz_y\nYour entry (layers): @write_layers";


		# DATASETS
		$input_datasets = get_input( "Which datasets (1 = yes, 0 = no, separated by comma)? (choice: @datavectorname)", $input_datasets, qr/[0-9, ]/ );
		my @write_datasets = map { 0 + $_ ? ($_/$_) : 0 + $_ } ( split /,/, $input_datasets )[0..3];
		$write_datasets[$_] = 0 for( $max_dataset + 1..2 );
		my $i = 0; print STDERR "Your entry: @write_datasets = " . join ' ', grep { $write_datasets[$i++] } @datavectorname;


		# FILENAME
		$input_filename = get_input( "Enter filename (\"%d\": diamonds; \"%l\": first layer; \"%i\" consecutive numbers)?", $input_filename, qr/[^[:space:]]/ );
		$filename_template = $input_filename;


		# ASSEMBLE INFORMATION
		@outputfiles = ( { 'diamonds_output' => [@write_diamonds], 'layers_output' => [@write_layers], 'data_output' => [@write_datasets], 'horiz_x_output' => [@write_horiz_x], 'horiz_y_output' => [@write_horiz_y] } );

		
		# WRITE OUTPUT FILE
		&{{ p => \&write_paraview_binary, g => \&write_gocad, x => \&write_xyz }->{$input_fileformat}};

		$output_filenumber++;
	}

}


###########################################################################
# EXPORT ROUTINES
###########################################################################

sub write_paraview_binary {
	my $fh;

	for my $outputfile (@outputfiles)
	{
		(my $output_filename = $filename_template) =~ s/%([dli])/{d => join( '', $outputfile->{'diamonds_output'}[0] ), l => $outputfile->{'layers_output'}[0], i => $output_filenumber }->{$1}/ge;
		open( $fh, '>', "$output_filename.vtk" );
		print STDERR "Writing to filename: $output_filename.vtk";
			
		my $numberofgridpoints = @{$outputfile->{'horiz_x_output'}} * @{$outputfile->{'horiz_y_output'}} * @{$outputfile->{'layers_output'}} * @{$outputfile->{'diamonds_output'}};

		# open Paraview output, write header
		print $fh "# vtk DataFile Version 3.0";
		print $fh "TERRA output";
		print $fh "binary";
		print $fh "DATASET UNSTRUCTURED_GRID";
		print $fh "POINTS $numberofgridpoints float";

		print STDERR "writing grid points ...";

		# write paraview grid
		for my $diamond ( @{$outputfile->{'diamonds_output'}} )
		{
			for my $layer ( @{$outputfile->{'layers_output'}} )
			{
				for my $i (@{$outputfile->{'horiz_x_output'}})
				{
					for my $j (@{$outputfile->{'horiz_y_output'}})
					{
						# scale grid coordinates with depth (given in meters), normalized to 1
						printf $fh '%s', pack( 'f>*', map {  $_ * ($layers[$layer] / 6.37E6) } @{$xg[$diamond][$i][$j]} );
					}
				}
			}
		}

		print $fh '';	# additional line break for binary

		print STDERR "writing connectivity ...";

		# write connectivity
		print $fh "CELLS " . (@{$outputfile->{'diamonds_output'}} * (@{$outputfile->{'layers_output'}} - 1) * (@{$outputfile->{'horiz_x_output'}} - 1) * (@{$outputfile->{'horiz_y_output'}} - 1) * 2) . " " . (7 * (@{$outputfile->{'diamonds_output'}} * (@{$outputfile->{'layers_output'}} - 1) * (@{$outputfile->{'horiz_x_output'}} - 1) * (@{$outputfile->{'horiz_y_output'}} - 1) * 2));

		for my $diamond (0..@{$outputfile->{'diamonds_output'}} - 1 )
		{
			for my $layer (0..@{$outputfile->{'layers_output'}} - 2 )
			{
				for my $i (0..@{$outputfile->{'horiz_x_output'}} - 2 )
				{
					for my $j (0..@{$outputfile->{'horiz_y_output'}} - 2 )
					{
						my $offset = $diamond * @{$outputfile->{'layers_output'}} * @{$outputfile->{'horiz_x_output'}} * @{$outputfile->{'horiz_y_output'}} + $layer * @{$outputfile->{'horiz_x_output'}} * @{$outputfile->{'horiz_y_output'}} + $i * @{$outputfile->{'horiz_y_output'}} + $j;
						my $offset_neighbor                     = $offset + 1;
						my $offset_next_row                     = $offset + @{$outputfile->{'horiz_y_output'}};
						my $offset_neighbor_next_row            = $offset_next_row + 1;
						my $offset_next_layer                   = $offset + @{$outputfile->{'horiz_x_output'}} * @{$outputfile->{'horiz_y_output'}};
						my $offset_neighbor_next_layer          = $offset_next_layer + 1;
						my $offset_next_row_next_layer          = $offset_next_layer + @{$outputfile->{'horiz_y_output'}};
						my $offset_neighbor_next_row_next_layer = $offset_next_row_next_layer + 1;

						printf $fh '%s', pack( 'l>*', 6, $offset, $offset_neighbor, $offset_next_row, $offset_next_layer, $offset_neighbor_next_layer, $offset_next_row_next_layer, 6, $offset_neighbor, $offset_neighbor_next_row, $offset_next_row, $offset_neighbor_next_layer, $offset_neighbor_next_row_next_layer, $offset_next_row_next_layer );
					}
				}
			}
		}
		
		print $fh '';	# additional line break for binary

		# all cells are trigonal prisms => cell type 13
		print $fh "CELL_TYPES " . (@{$outputfile->{'diamonds_output'}} * (@{$outputfile->{'layers_output'}} - 1) * (@{$outputfile->{'horiz_x_output'}} - 1) * (@{$outputfile->{'horiz_y_output'}} - 1) * 2);
		printf $fh '%s', pack( 'l>', 13 ) for( 1..(@{$outputfile->{'diamonds_output'}} * (@{$outputfile->{'layers_output'}} - 1) * (@{$outputfile->{'horiz_x_output'}} - 1) * (@{$outputfile->{'horiz_y_output'}} - 1) * 2) );

		print $fh "POINT_DATA $numberofgridpoints";
		
		print $fh '';	# additional line break for binary

		for my $dataset (0..3)
		{
			next unless( $outputfile->{'data_output'}[$dataset] );
			print STDERR "writing data $dataset ($datavectorname[$dataset]) ...";

			print $fh (3 == $datavectorsize[$dataset] ? "VECTORS $datavectorname[$dataset] float" : "SCALARS $datavectorname[$dataset] float $datavectorsize[$dataset]\nLOOKUP_TABLE default");

			for my $diamond ( @{$outputfile->{'diamonds_output'}} )
			{
				for my $layer ( @{$outputfile->{'layers_output'}} )
				{
					for my $i (@{$outputfile->{'horiz_x_output'}})
					{
						for my $j (@{$outputfile->{'horiz_y_output'}})
						{
							if( 3 == $dataset )
							{
								printf( $fh '%s', pack( 'f>', $data[ get_index( 0, $layer, $diamond, $i, $j, 0 ) ] - $mean[$layer] ) );
							}
							else
							{
								printf( $fh '%s', pack( 'f>*', map { $data[ get_index( $dataset, $layer, $diamond, $i, $j, $_ ) ] } ( 0..$datavectorsize[$dataset] - 1 ) ) );
							}
						}
					}
				}
			}
			
			print $fh '';	# additional line break for binary

		}

		close( $fh );
	}
}
	
sub write_gocad {
	my $fh;

	for my $outputfile (@outputfiles)
	{
		(my $output_filename = $filename_template) =~ s/%([dli])/{d => join( '', $outputfile->{'diamonds_output'}[0] ), l => $outputfile->{'layers_output'}[0], i => $output_filenumber }->{$1}/ge;
		open( $fh, '>', "$output_filename.ts" );
		print STDERR "Writing to filename: $output_filename.ts";
				
		my $numberofgridpoints = @{$outputfile->{'horiz_x_output'}} * @{$outputfile->{'horiz_y_output'}} * @{$outputfile->{'layers_output'}} * @{$outputfile->{'diamonds_output'}};

		# open Gocad output, write header
		print $fh "GOCAD TSolid 1 ";
		print $fh "HEADER {";
		print $fh "name:TERRA";
		print $fh "}";
		$l = 0; print $fh "PROPERTIES " . join ' ', grep { $outputfile->{'data_output'}[$l++] == 1 } @datavectorname;
		$l = 0; print $fh "ESIZES " . join ' ', grep { $outputfile->{'data_output'}[$l++] == 1 } @datavectorsize;
		$l = 0; print $fh "UNITS " . join ' ', grep { $outputfile->{'data_output'}[$l++] == 1 } ('K', 'm/s', 'Pa', 'K');
		print $fh "GOCAD_ORIGINAL_COORDINATE_SYSTEM";
		print $fh "NAME Default";
		print $fh "AXIS_NAME \"X\" \"Y\" \"Z\"";
		print $fh "AXIS_UNIT \"m\" \"m\" \"m\"";
		print $fh "ZPOSITIVE Elevation";
		print $fh "END_ORIGINAL_COORDINATE_SYSTEM";
		print $fh "TVOLUME";

		print STDERR "writing grid points ...";

		$k = 0;
		# write gocad grid
		for my $diamond ( @{$outputfile->{'diamonds_output'}} )
		{
			for my $layer ( @{$outputfile->{'layers_output'}} )
			{
				for my $i (@{$outputfile->{'horiz_x_output'}})
				{
					for my $j (@{$outputfile->{'horiz_y_output'}})
					{
						$l = 0;

						# scale grid coordinates with depth (given in meters), normalized to 1
						print $fh join " ", ("PVRTX", ++$k, map { $_ * ($layers[$layer] / 6.37E6) } @{$xg[$diamond][$i][$j]}, map { @$_ } grep { $outputfile->{'data_output'}[$l++] == 1 } ( [ map { $data[ get_index( 0, $layer, $diamond, $i, $j, $_ ) ] } ( 0..$datavectorsize[0] - 1 ) ], [ map { $data[ get_index( 1, $layer, $diamond, $i, $j, $_ ) ] } ( 0..$datavectorsize[1] - 1 ) ], [ map { $data[ get_index( 2, $layer, $diamond, $i, $j, $_ ) ] } ( 0..$datavectorsize[2] - 1 ) ], [ $data[ get_index( $dataset, $layer, $diamond, $i, $j, 0 ) ] - $mean[$layer] ] ) ); 
					}
				}
			}
		}

		print STDERR "writing connectivity ...";

		# write connectivity
		for my $diamond ( @{$outputfile->{'diamonds_output'}} )
		{
			 for my $layer (0..@{$outputfile->{'layers_output'}} - 2 )
			{
				for my $i (0..@{$outputfile->{'horiz_x_output'}} - 2 )
				{
					for my $j (0..@{$outputfile->{'horiz_y_output'}} - 2 )
					{
						local $, = " ";

						my $offset = $diamond * @{$outputfile->{'layers_output'}} * @{$outputfile->{'horiz_x_output'}} * @{$outputfile->{'horiz_y_output'}} + $layer * @{$outputfile->{'horiz_x_output'}} * @{$outputfile->{'horiz_y_output'}} + $i * @{$outputfile->{'horiz_y_output'}} + $j + 1;
						my $offset_neighbor                     = $offset + 1;
						my $offset_next_row                     = $offset + @{$outputfile->{'horiz_y_output'}};
						my $offset_neighbor_next_row            = $offset_next_row + 1;
						my $offset_next_layer                   = $offset + @{$outputfile->{'horiz_x_output'}} * @{$outputfile->{'horiz_y_output'}};
						my $offset_neighbor_next_layer          = $offset_next_layer + 1;
						my $offset_next_row_next_layer          = $offset_next_layer + @{$outputfile->{'horiz_y_output'}};
						my $offset_neighbor_next_row_next_layer = $offset_next_row_next_layer + 1;

						# first prism:
						# upper triangle and point in next layer (below first point of triangle)
						print $fh "TETRA", $offset, $offset_neighbor, $offset_next_row, $offset_next_layer;

						# two points in each layer
						print $fh "TETRA", $offset_neighbor, $offset_next_row, $offset_next_layer, $offset_neighbor_next_layer;

						# lower triangle and point in upper layer (in next row)
						print $fh "TETRA", $offset_next_row, $offset_next_layer, $offset_neighbor_next_layer, $offset_next_row_next_layer;


						# second prism:
						# upper triangle and point in next layer (below first point of triangle)
						print $fh "TETRA", $offset_neighbor, $offset_neighbor_next_row, $offset_next_row, $offset_neighbor_next_layer;

						# two points in each layer
						print $fh "TETRA", $offset_neighbor_next_row, $offset_next_row, $offset_neighbor_next_layer, $offset_neighbor_next_row_next_layer;

						# lower triangle and point in upper layer (above second point of triangle)
						print $fh "TETRA", $offset_next_row, $offset_neighbor_next_layer, $offset_neighbor_next_row_next_layer, $offset_next_row_next_layer;
					}
				}
			}
		}

		print $fh "END";
		
		close( $fh );
	}
}


sub write_xyz {
	my $fh;

	for my $outputfile (@outputfiles)
	{
		(my $output_filename = $filename_template) =~ s/%([dli])/{d => join( '', $outputfile->{'diamonds_output'}[0] ), l => $outputfile->{'layers_output'}[0], i => $output_filenumber }->{$1}/ge;
		open( $fh, '>', "$output_filename.xyz" );
		print STDERR "Writing to filename: $output_filename.xyz";

		my $numberofgridpoints = @{$outputfile->{'horiz_x_output'}} * @{$outputfile->{'horiz_y_output'}} * @{$outputfile->{'layers_output'}} * @{$outputfile->{'diamonds_output'}};

		print STDERR "writing grid points ...";

		for my $diamond ( @{$outputfile->{'diamonds_output'}} )
		{
			for my $layer ( @{$outputfile->{'layers_output'}} )
			{
				my $depth = 6.37E6 - $layers[$layer];

				for my $i (@{$outputfile->{'horiz_x_output'}})
				{
					for my $j (@{$outputfile->{'horiz_y_output'}})
					{
						$l = 0;

						my $lon = rad2deg( atan2( $xg[$diamond][$i][$j][1], $xg[$diamond][$i][$j][0] ) );
						my $length = sqrt( $xg[$diamond][$i][$j][0]**2 + $xg[$diamond][$i][$j][1]**2 );
						my $lat = rad2deg( atan2( ( $length ), ( $xg[$diamond][$i][$j][2] ) ) );
						$lat += 180 if( $lat < 0 );
						$lat = 90 - $lat;

						print $fh join "\t", ($lon, $lat, $depth, map { @$_ } grep { $outputfile->{'data_output'}[$l++] == 1 } ( [ map { $data[ get_index( 0, $layer, $diamond, $i, $j, $_ ) ] } ( 0..$datavectorsize[0] - 1 ) ], [ map { $data[ get_index( 1, $layer, $diamond, $i, $j, $_ ) ] } ( 0..$datavectorsize[1] - 1 ) ], [ map { $data[ get_index( 2, $layer, $diamond, $i, $j, $_ ) ] } ( 0..$datavectorsize[2] - 1 ) ], [ $data[ get_index( $dataset, $layer, $diamond, $i, $j, 0 ) ] - $mean[$layer] ] ) ); 
					}
				}
			}
		}

		close( $fh );
	}
}

__DATA__
1	# number of datasets to be read
#
# First file
#
p	# type = Paraview
0	# diamond = 1
0-64	# x direction
0-64	# y direction
0-32	# layers
2	# downsampling factor
1,0,0,1	# datasets
out%i	# filename without extension
#
# Second file
#
p
1
0-64
0-64
0-32
2
1,0,0,1
out%i
q	# quit
__END__
==pod
=head1 NAME

terra_export.pl - An export program for TERRA output files

=head1 SYNOPSIS

B<terra_export.pl> inputfiles ...

B<terra_export.pl> B<-e>

B<terra_export.pl> B<-h>

=head1 DESCRIPTION

terra_export.pl is an export program which can export TERRA simulations into several output formats for visualization:

=over

=item * Paraview (binary *.vtk unstructured grid format)

=item * Gocad (TSolid *.ts format)

=item * xyz format (lon, lat, data values, separated by tabulator)

=back

The program deduces the grid size from the number of input files and the header of the first input file. Thus, always the whole model has to be loaded. (If you are really just interested in a part of the dataset, specify the right number of empty files on the command line such that the total number of files is correct and that at least the first file contains a correct header.) The program generates the grid and connectivity and asks the user interactively for output options.

However, the program can also be run in a scripted mode. For this, the script is piped into STDIN. The script contains the same keypresses as in the interactive version. Comments (starting with the '#' sign) are ignored.

An example script can be requested by using the program with the parameter I<-h>.

Scripts can also be recorded by piping STDOUT into a file. (All other output is written to STDERR.)

=head1 OPTIONS

B<-e>
	print the example script and exit the program

B<-h>
	show help (this page)

=head1 EXAMPLES

B<terra_export.pl> F<c300.00*.01>

=over

run interactively

=back

B<terra_export.pl> F<c300.00*.01.gz>

=over

run interactively, read compressed input files; for this, I<zcat> must be available

=back

B<terra_export.pl> -h

=over

show help (this page)

=back

B<terra_export.pl> -e
	
=over

show example script

=back

B<terra_export.pl> -e > F<examplescript.txt>

=over

save example script

=back

B<terra_export.pl> F<c300.00*.01> > F<script.txt>

=over

run interactively, record inputs to script.txt

=back

B<terra_export.pl> F<c300.00*.01> < F<script.txt>

=over

run non-interactively, controlled by script

=back

B<terra_export.pl> -e | B<terra_export.pl> F<c300.00*.01>

=over

apply example script to dataset

=back

B<terra_export.pl> F<c300.0000.01> $(seq 1 31)

=over

read only first partition of an MT64 simulation, replace rest of input files by dummy files (numbers from 1 to 31)

=back

B<perl> -pe '$.-1or$l=int($_/10+9);$.-$l||exit' F<c300.0016.03> > F<header.txt>; B<terra_export.pl> F<header.txt> $(seq 1 15) F<c300.0016.03> $(seq 17 31)

=over

create dummy header as first file; read only 17th partition of an MT64 simulation, replace rest of input files by dummy files

=back

=head1 BUGS

The input syntax of the layers has still some bugs.

=head1 AUTHOR

Christoph Moder <moder@geophysik.uni-muenchen.de>

=head1 HISTORY

August 2011: first version released to http://earthmodels.org/

October 2011: new features: downsampling, variable horizontal size, scripting

This version: 2011-10-20

=cut