#!/usr/bin/perl -wanl

# USAGE:
# $0 polygons.xy [nt] > cfile [3> paraviewfile.vtk]
#
# ALL TIME STEPS (16 parallel processes):
# { echo -n "Run started: "; date; for((i=0; i<=200; i++)); do echo -n "Process $i: "; date; ( ./findplate_xy_gplates.pl Polygons.platepolygons.$i.platepolygons.$i.xy > platemap_$i.dat 2> /dev/null & ) ; while( [ $(ps ux | grep '[f]indplate_xy_gplates' | wc -l) -ge 16 ] ); do sleep 5; done; done; while( [ $(ps ux | grep '[f]indplate_xy_gplates' | wc -l) -gt 0 ] ); do sleep 5; done; echo -n "Run ended: "; date; } | tee findplate.log
#
# AUTHOR OF THIS SCRIPT: Christoph Moder <moder@geophysik.uni-muenchen.de>
# ACTUAL ALGORITHM: taken from GPlates 0.9.9.1 (http://gplates.org/), file src/maths/PointInPolygon.cc

use constant NT => ( $ARGV[-1] =~ /^[2468]|[1-9][0-9]*[02468]$/ ? pop( @ARGV ) : 256 );	# number of grid subdivisions; default value = 256, or specified as second parameter on command line
use constant NR => NT / 2;	# number of radial layers, usually NT/2
use constant PI => atan2(0, -1);
use constant NORTH => 0;
use constant SOUTH => 10000;

sub is_point_in_on_out_counter {
	my @point = @_[0..1];
	my $count_north = \$_[2];
	my $count_south = \$_[3];
	my @polygon = @_[4..$#_];
	my $polypoints = scalar( @polygon );

	# local tmp vars
	# using Real lets us use operator== 
	my ($W, $E, $S, $N, $x_lat, $dlon, $lon, $lon1, $lon2);

	# re-set number of crossings
	$$count_south = 0;
	$$count_north = 0;

	# Compute meridian through P and count all the crossings with 
	# segments of polygon boundary 

	for( my $i = 0; $i <= $#polygon; $i++ )
	{
		@v1 = @{$polygon[$i]};
		@v2 = @{$polygon[($i + 1) % $polypoints]};

		# Copy the two vertex longitudes 
		# since we need to mess with them
		$lon1 = $v1[0];
		$lon2 = $v2[0];

		# delta in lon
		$dlon = $lon2 - $lon1;

		# Jumped across Greenwhich going westward or eastward
		( ( $dlon > 180.0 ) and $lon2 -= 360.0 ) or ( ( $dlon < -180.0 ) and $lon1 -= 360.0 );

		# set lon limits for this segment
		if ($lon1 <= $lon2) 
		{	
			# segment goes W to E (or N-S)
			$W = $lon1;
			$E = $lon2;
		}
		else 
		{			
			# segment goes E to W
			$W = $lon2;
			$E = $lon1;
		}
		
		# Local copy of plon, adjusted given the segment lon range
		$lon = $point[0];	

		# Make sure we rewind way west for starters
		$lon -= 360.0 while ($lon > $W);

		# Then make sure we wind to inside the lon range or way east
		$lon += 360.0 while ($lon < $W);

		# Not crossing this segment
		next if ($lon > $E); # to next vertex

		# Special case of N-S segment: does P lie on it?
		if ($dlon == 0.0) 
		{	
			if ( $v2[1] < $v1[1] ) 
			{	
				# Get N and S limits for segment
				$S = $v2[1];
				$N = $v1[1];
			}
			else
			{
				$N = $v2[1];
				$S = $v1[1];
			}

			# P is not on this segment
			next if ( $point[1] < $S || $point[1] > $N ); # to next vertex

			# P is on segment boundary; we are done
			return (1);
		}

		# Calculate latitude at intersection
		$x_lat = $v1[1] + ( ($v2[1] - $v1[1]) / ($lon2 - $lon1) ) * ($lon - $lon1);

		# P is on S boundary
		return 1 if ( $x_lat == $point[1] );

		# Only allow cutting a vertex at end of a segment 
		# to avoid duplicates 
		next if ($lon == $lon1);

		# Is cut is north or south of P?
		($x_lat > $point[1]) and ++$$count_north or ++$$count_south;

	} # end of loop over vertices

	return (0);	
}

sub rad2deg {
	return ($_[0] * 180 / PI);
}

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

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 );

	for( my $diamond = 0; $diamond < $nd; $diamond++ )
	{
		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; $k < $lvt; $k++ )
		{
			my $m = 2**$k;
			my $l = NT / $m;

			# rows of diamond
			for( my $j1 = 0; $j1 <= $m; $j1++ )
			{
				for( my $j2 = 0; $j2 < $m; $j2++ )
				{
					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][$_] for( 0..2 );

					# 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; $j1 <= $m; $j1++ )
			{
				for( my $j2 = 0; $j2 < $m; $j2++ )
				{
					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][$_] for( 0..2 );

					# 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; $j1 < $m; $j1++ )
			{
				for( my $j2 = 0; $j2 < $m; $j2++ )
				{
					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][$_] for( 0..2 );
					
					# 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 );
				}
			}
		}
	}
}

# adjusts the minimum or maximum value if the given value is below or above
sub minmax {
	( ( $_[2] < $_[0] ) and ( $_[0] = $_[2] ) ) or ( ( $_[2] > $_[1] ) and ( $_[1] = $_[2] ) );
}

# read polygons
if(/^> /) {
	# close previous polygons
	if( defined( $current_plate ) )
	{
		for( (NORTH, SOUTH) )
		{
			if( $polygons{$current_plate + $_}->[-1] )
			{
				$dlon = $polygons{$current_plate + $_}->[-1]->[0] - $polygons{$current_plate + $_}->[0]->[0];
				$dlon = abs( 360.0 - abs( $dlon ) ) * (-$dlon <=> 0) if( abs( $dlon ) > 180.0 );
				$bbox{$current_plate + $_}->{'lonsum'} += $dlon;
			}
		}
	}
	
	# new plate: extract id
	$current_plate = (/;.+?_(\d{3})_ /)[0];	# behind first semicolon: three digits between underscores, then space

	for( (NORTH, SOUTH) )
	{
		$bbox{$current_plate + $_}->{'minlat'} = 90;
		$bbox{$current_plate + $_}->{'maxlat'} = -90;
		$bbox{$current_plate + $_}->{'minlon'} = 180;
		$bbox{$current_plate + $_}->{'maxlon'} = -180;
		$bbox{$current_plate + $_}->{'pole'} = 0;
		$bbox{$current_plate + $_}->{'lonsum'} = 0;
	}
	undef $prev_lon, $prev_lat if( defined( $prev_lon ) );
}
# coordinates: store them in array, determine min/max
elsif( $#F == 1 )
{
	my $new_hemisphere = ($F[1] < 0) * SOUTH;

	# we know that first and last point in the polygon is identical
	if( defined( $prev_lon ) )
	{
		my $old_hemisphere = ($prev_lat < 0) * SOUTH;

		$dlon = $prev_lon - $F[0];
		$dlon = abs( 360.0 - abs( $dlon ) ) * (-$dlon <=> 0) if( abs( $dlon ) > 180.0 );

		# crossed the equator?
		if( $old_hemisphere != $new_hemisphere )
		{
			# determine longitude of intersection
			$lon_middle = $prev_lon + $dlon / ($prev_lat - $F[1]) * (-$prev_lat);
			$lon_middle = abs( 360.0 - abs( $lon_middle ) ) * (-$lon_middle <=> 0) if( abs( $lon_middle ) > 180.0 );

			# add point of equator intersection to previous plate
			push( @{ $polygons{$current_plate + $old_hemisphere} }, [ $lon_middle, 0 ] );
			minmax( $bbox{$current_plate + $old_hemisphere}->{'minlat'}, $bbox{$current_plate + $old_hemisphere}->{'maxlat'}, 0 );
			minmax( $bbox{$current_plate + $old_hemisphere}->{'minlon'}, $bbox{$current_plate + $old_hemisphere}->{'maxlon'}, $lon_middle );

			$dlon = $prev_lon - $lon_middle;
			$dlon = abs( 360.0 - abs( $dlon ) ) * (-$dlon <=> 0) if( abs( $dlon ) > 180.0 );
			$bbox{$current_plate + $old_hemisphere}->{'lonsum'} += $dlon;

			# add longitude difference of equator segment of following polygon, if polygon of this hemisphere exists
			if( $polygons{$current_plate + $new_hemisphere}->[-1] )
			{
				$prev_lon_thisplate = $polygons{$current_plate + $new_hemisphere}->[-1]->[0];
				$dlon = $prev_lon_thisplate - $lon_middle;
				$dlon = abs( 360.0 - abs( $dlon ) ) * (-$dlon <=> 0) if( abs( $dlon ) > 180.0 );
				$bbox{$current_plate + $new_hemisphere}->{'lonsum'} += $dlon;
			}

			# add point of equator intersection to following plate
			push( @{ $polygons{$current_plate + $new_hemisphere} }, [ $lon_middle, 0 ] );
			minmax( $bbox{$current_plate + $new_hemisphere}->{'minlat'}, $bbox{$current_plate + $new_hemisphere}->{'maxlat'}, 0 );
			minmax( $bbox{$current_plate + $new_hemisphere}->{'minlon'}, $bbox{$current_plate + $new_hemisphere}->{'maxlon'}, $lon_middle );

			$dlon = $lon_middle - $F[0];
			$dlon = abs( 360.0 - abs( $dlon ) ) * (-$dlon <=> 0) if( abs( $dlon ) > 180.0 );
			$bbox{$current_plate + $new_hemisphere}->{'lonsum'} += $dlon;

			if( not defined( $bbox{$current_plate + $new_hemisphere}->{'firstpoint'} ) )
			{
				$bbox{$current_plate + $new_hemisphere}->{'firstpoint'} = [ $lon_middle, 0 ];
			}

		}
		else
		{
			$bbox{$current_plate + $new_hemisphere}->{'lonsum'} += $dlon;
		}
	}
	
	push( @{ $polygons{$current_plate + $new_hemisphere} }, [ $F[0], $F[1] ] );
	
	minmax( $bbox{$current_plate + $new_hemisphere}->{'minlat'}, $bbox{$current_plate + $new_hemisphere}->{'maxlat'}, $F[1] );
	minmax( $bbox{$current_plate + $new_hemisphere}->{'minlon'}, $bbox{$current_plate + $new_hemisphere}->{'maxlon'}, $F[0] );

	$bbox{$current_plate + $new_hemisphere}->{'firstpoint'} = [ $F[0], $F[1] ] unless( defined( $bbox{$current_plate + $new_hemisphere}->{'firstpoint'} ) );

	$prev_lon = $F[0];
	$prev_lat = $F[1];
}

END {
	$point = 0;
	$plate = 0;
	$nd = 10;	# number of diamonds
	
	print STDERR "polygons read";
	$numberofgridpoints = (NT + 1)**2 * $nd;
	
	# dermine if the platepolygon contains the pole
	( abs( abs( $bbox{$_}->{'lonsum'} ) - 360 ) < 1.0e-8 ) and $bbox{$_}->{'pole'} = ( ( abs( $bbox{$_}->{'maxlat'} ) > abs( $bbox{$_}->{'minlat'} ) ) ? ( $bbox{$_}->{'maxlat'} <=> 0 ) : ( $bbox{$_}->{'minlat'} <=> 0 ) ) for ( keys %polygons );
	
	# write TERRA header
	# first: grid resolution (radial layers, diamond subdivisions)
	printf( "%5d%5d\n", NR, NT );

	# four lines of comment, for the case name
	print "comment 1";
	print "comment 2";
	print "comment 3";
	print "comment 4";
	
	# radial layer depth
	printf( "%15.8E%s", 0, ($_ % 10 ? "" : "\n") ) for( 1..(NT / 2 + 1) );

	print "";

	# property array
	printf( "%15.8E%s", 0, ($_ % 10 ? "" : "\n") ) for( 1..20 );
	
	# generate grid
	grdgen();
	print STDERR "grid generated";

	# open Paraview output, write header
	open( PARAVIEW, ">&3" ) or open( PARAVIEW, ">", "/dev/null" );
	print PARAVIEW "# vtk DataFile Version 3.0";
	print PARAVIEW "TERRA plate map";
	print PARAVIEW "ASCII";
	print PARAVIEW "DATASET POLYDATA";
	print PARAVIEW "POINTS $numberofgridpoints float";
	
	# write paraview grid
	for( my $diamond = 0; $diamond < $nd; $diamond++ )
	{
		for( my $i = 0; $i <= NT; $i++ )
		{
			for( my $j = 0; $j <= NT; $j++ )
			{
				print PARAVIEW "@{$xg[$diamond][$i][$j]}";
			}
		}
	}

	# write connectivity
	print PARAVIEW "POLYGONS " . (10 * NT**2 * 2) . " " . (4 * (10 * NT**2 * 2));
	for( my $diamond = 0; $diamond < $nd; $diamond++ )
	{
		for( my $i = 0; $i < NT; $i++ )
		{
			for( my $j = 0; $j < NT; $j++ )
			{
				local $, = " ";
				my $offset = $diamond * (NT + 1)**2 + $i * (NT + 1) + $j;
				print PARAVIEW 3, $offset, $offset + 1, $offset + (NT + 1);
				print PARAVIEW 3, $offset + 1, $offset + (NT + 1) + 1, $offset + (NT + 1);
			}
		}
	}
 	print PARAVIEW "POINT_DATA $numberofgridpoints";
	print PARAVIEW "SCALARS plate_id float";
	print PARAVIEW "LOOKUP_TABLE default";

	# do the actual work: assign plate ids to grid points
	for( my $diamond = 0; $diamond < $nd; $diamond++ )
	{
		for( my $i = 0; $i <= NT; $i++ )
		{
			for( my $j = 0; $j <= NT; $j++ )
			{
				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;
		
				$prev_plate = $plate;
				$plate = 0;	# not found
				
				# test all plates; first: the one of the previous grid point; then: plates with few polygon points first, so the expensive tests with many polygon points can be avoided if possible
				for $current_plate ( sort { $a == $prev_plate ? -1 : $b == $prev_plate ? +1 : scalar(@{$polygons{$a}}) <=> scalar(@{$polygons{$b}}) } keys %polygons )
				{
					# Algorithm:
					#
					# Case 1: The polygon S contains a geographical pole
					#	   a) if P is beyond the far latitude then P is outside
					#	   b) Compute meridian through P and count intersections:
					#		odd: P is outside; even: P is inside
					#
					# Case 2: S does not contain a pole
					#	   a) If P is outside range of latitudes then P is outside
					#	   c) Compute meridian through P and count intersections:
					#		odd: P is inside; even: P is outside
					#
					# In all cases, we check if P is on the outline of S
					
					#counters for the number of crossings of a meridian through p and and the segments of this polygon
					my $count_north = 0;
					my $count_south = 0;

					if( $bbox{$current_plate}->{'pole'} ) 
					{	
						# Case 1 of an enclosed polar cap

						# N polar cap
						if( $bbox{$current_plate}->{'pole'} == 1 )
						{	
							# South of a N polar cap
							next if( $lat < $bbox{$current_plate}->{'minlat'} );

							# Clearly inside of a N polar cap
							$plate = $current_plate and last if( $lat > $bbox{$current_plate}->{'maxlat'} );
						}

						# S polar cap
						if( $bbox{$current_plate}->{'pole'} == -1 )
						{	
							# North of a S polar cap
							next if( $lat > $bbox{$current_plate}->{'maxlat'} );

							# North of a S polar cap
							$plate = $current_plate and last if( $lat < $bbox{$current_plate}->{'minlat'} );
						}
					
						# Tally up number of intersections between polygon 
						# and meridian through p 
						
						$plate = $current_plate and last if( is_point_in_on_out_counter( $lon, $lat, $count_north, $count_south, @{$polygons{$current_plate}}) );	# Found P is on S; return (POINT_ON_POLYGON);	
					
						$plate = $current_plate and last if( $bbox{$current_plate}->{'pole'} == 1 and $count_north % 2 == 0);

						$plate = $current_plate and last if( $bbox{$current_plate}->{'pole'} == -1 and $count_south % 2 == 0);
					
						next;
					}
					
					# Here is Case 2.

					# First check latitude range 
					next if( $lat < $bbox{$current_plate}->{'minlat'} or $lat > $bbox{$current_plate}->{'maxlat'} );
					
					# Longitudes are tricker and are tested with the tallying of intersections 
					
					$plate = $current_plate and last if( is_point_in_on_out_counter( $lon, $lat, $count_north, $count_south, @{$polygons{$current_plate}}) );	# Found P is on S; return (POINT_ON_POLYGON);	
					$plate = $current_plate and last if( $count_north % 2 );	# point inside polygon
					
					# Nothing triggered the tests; we are outside
				}

				$could_use_prev++ if( $plate == $prev_plate );
				
				$plate = $prev_plate unless( $plate );	# use previous plate id if none has been found

				printf( "%10.3f%s", ($plate < SOUTH ? $plate : $plate - SOUTH), (++$point % 15 ? "" : "\n") );
				print PARAVIEW ($plate < SOUTH ? $plate : $plate - SOUTH);

				# print progress in percent on STDERR
				print STDERR int($point / $numberofgridpoints * 100 + 0.5) . "%" unless( $point % ($numberofgridpoints / 100) );
			}
		}
	}

	print STDERR "efficiency: $could_use_prev of $numberofgridpoints = " . int($could_use_prev * 100 / $numberofgridpoints + 0.5) . "%";
}