#!/usr/bin/perl -anl

use constant PI => atan2(0, -1);

sub deg2rad {  $_[0] / 180 * PI }

BEGIN {
	# connect last longitude with first?
	$is_sphere = 1;

	# convert lon/lat to cartesian (= true) or assume that the coordinates are already cartesian?
	$polar2cart = 1;

	# should there be an exaggeration in radial direction (or z direction if non-spherical)?
	$exaggerate = 0;
	$radius = 6370E3;

	$, = "\t";
	$" = "\t";

	print STDERR "reading data ...";
}

($lon, $lat, $altitude) = map { 0 + $_ } @F;

$xwidth{$lon}++;
$ywidth{$lat}++;
$point[$.-1] = [ @F ];

END {
	$gridpoints = $.;
	print STDERR "writing grid points ($gridpoints points, " . keys( %xwidth ) . "x" . keys( %ywidth ) . ") ...";
	print STDERR "x direction reversed" if( $point[0]->[0] > $point[-1]->[0] );
	print STDERR "y direction reversed" if( $point[0]->[1] > $point[-1]->[1] );

	# output vtk header
	print "# vtk DataFile Version 3.0";
	print "xyz2vtkpoly.pl";
	print "ASCII";
	print "DATASET POLYDATA";
	print "POINTS $gridpoints float";

	# output grid points
	if( $polar2cart )
	{
		print "@{$_}" for( map { my $r = ($radius + $_->[2] * $exaggerate) / $radius; my $x = $r * cos( deg2rad($_->[0]) ) * cos( deg2rad($_->[1]) ); my $y = $r * sin( deg2rad($_->[0]) ) * cos( deg2rad($_->[1]) ); my $z = $r * sin( deg2rad($_->[1]) ); [$x, $y, $z]; } @point );
	}
	else
	{
		print "@{$_}" for( map { [$_->[0], $_->[1], $_->[2] * $exaggerate] } @point );
	}

	# output polygons
	print STDERR "writing connectivity ...";
	print "POLYGONS " . ($is_sphere ? ($gridpoints - keys( %xwidth )) . " " . 5 * ($gridpoints - keys( %xwidth ) ) : ($gridpoints - keys( %xwidth ) - keys( %ywidth ) + 1) . " " . 5 * ($gridpoints - keys( %xwidth ) - keys( %ywidth ) + 1) ); 

	for my $lat ( ($point[0]->[1] < $point[-1]->[1]) ? 0..keys( %ywidth ) - 2 : reverse 0..keys( %ywidth ) - 2 )
	{
		for my $lon ( ($point[0]->[0] < $point[-1]->[0]) ? 0..keys( %xwidth ) - 2 : reverse 0..keys( %xwidth ) - 2 )
		{
			print 4, ($lat * keys( %xwidth ) + $lon), ( ($lat + 1) * keys( %xwidth ) + $lon), ( ($lat + 1) * keys( %xwidth ) + $lon + 1), ($lat * keys( %xwidth ) + $lon + 1);
		}
			
		$is_sphere and print 4, ($lat * keys( %xwidth ) + $lon),( ($lat + 1) * keys( %xwidth ) + $lon), ( ($lat + 1) * keys( %xwidth ) ), ($lat * keys( %xwidth ));
	}

	# output height data
	print STDERR "writing data ...";
	print "POINT_DATA $gridpoints";
	print "SCALARS altitude float";

	print for( map { $_->[2] } @point );
}