#!/usr/bin/perl

# 2010mar06a PMT : allow arbitary regular quadrilateral orientation
# 2006mar29a PMT : estimate area of a rectangle given shape model
#  (of Eros) and lat, lon of rectangle corners

use strict;
use warnings;
use POSIX qw /floor/;
use Math::Trig qw / :pi deg2rad rad2deg /;

sub fill_out_model
{
    my ($shape_p) = @_;

    for (my $lat = -90; $lat <= 90; $lat+=2) { # complete these first
	for (my $lon = 1; $lon < 360; $lon+=2) {
	    store($shape_p, $lat, $lon,
		  (retrieve($shape_p, $lat, $lon-1) +
		   retrieve($shape_p, $lat, $lon+1))/2.);
	}
    }
    for (my $lat = -89; $lat < 90; $lat+=2) { # across lat for these
	for (my $lon = 0; $lon <= 360; $lon++) {
	    store($shape_p, $lat, $lon,
		  (retrieve($shape_p, $lat-1, $lon) +
		   retrieve($shape_p, $lat+1, $lon))/2.);
	}
    }
}


sub integrate_old
{
    my ($s_p, $c_ap) = @_;
    # treat as aligned w/ lat/lon array
    warn "approximating as aligned in lat lon\n";
    my ($lat_min, $lat_max, $lon_min, $lon_max) = (91, -91, 361, -1);
    for (my $i = 0; $i < 4; $i++) {
	$lat_min = ${$c_ap->[$i]}->{lat} if (${$c_ap->[$i]}->{lat} < $lat_min);
	$lat_max = ${$c_ap->[$i]}->{lat} if (${$c_ap->[$i]}->{lat} > $lat_max);
	$lon_min = ${$c_ap->[$i]}->{lon} if (${$c_ap->[$i]}->{lon} < $lon_min);
	$lon_max = ${$c_ap->[$i]}->{lon} if (${$c_ap->[$i]}->{lon} > $lon_max);
    }
    print "integrating over latitudes $lat_min to $lat_max "
	."and longitudes $lon_min to $lon_max\n";
    my $dlat = ($lat_max - $lat_min) / 1000.;
    my $dlon = ($lon_max - $lon_min) / 1000.;
    my $area = 0.;
    for (my $lat = $lat_min; $lat <= $lat_max; $lat += $dlat) {
	my $dloncos = cos($lat * 3.14159 / 180.) * $dlon;
	for (my $lon = $lon_min; $lon <= $lon_max; $lon += $dlon) {
	    $area += interpolate($s_p, $lat, $lon)**2 * $dlat * $dloncos
		* 3.14159 / 180. * 3.14159 / 180.;
	}
    }
    return $area;
}


sub integrate
{
    my ( $shape_p, $corners_p ) = @_;
    die "only quadrilaterals allowed\n" unless scalar @$corners_p == 4;

    # enforce an order on the vertices.  CW from 0
    #  start by finding middle of quadrilateral
    my ( $clon, $clat ) = ( 0, 0 );
    for ( my $i = 0; $i < 4; $i++ ) {
	$clon += $corners_p->[$i]->{lon};
	$clat += $corners_p->[$i]->{lat};
    }
    $clon /= 4.0;
    $clat /= 4.0;

    #  find direction of each vertex from center
    for ( my $i = 0; $i < 4; $i++ ) {
	$corners_p->[$i]->{angle} = atan2(
	    $corners_p->[$i]->{lat} - $clat,
	    $corners_p->[$i]->{lon} - $clon );
	# force to be in range 0->pi2
	$corners_p->[$i]->{angle} += pi2 if $corners_p->[$i]->{angle} < 0;
	die "bad angle\n" if $corners_p->[$i]->{angle} < 0
	    or $corners_p->[$i]->{angle} > pi2;
    }
    # sort by that angle
    my @p =           sort { $a->{angle} <=> $b->{angle} } @$corners_p;
    my @lat_ordered = sort { $a->{lat}   <=> $b->{lat}   } @$corners_p;

    # left-most triangle is always lat_ordered[0->1]
    #  one of its sides is always that vector.  But the other side
    #  can be to either of the other points.
    my $iam;
    for ( my $i = 0; $i < 4; $i++ ) {
	$iam = $i if $p[$i] == $lat_ordered[0];  # _pointer_ comparison
    }
    die "cannot find myself\n" unless $iam;
    my $left = ($iam+3)%4;
    my $right = ($iam+5)%4;

    # slopes are 90d rotated because more effient to integrate over lon
    #  than over lat due to cos(lat) term
    my $slope_left = ( $p[$left]->{lon} - $p[$iam]->{lon} )
	/ ( $p[$left]->{lat} - $p[$iam]->{lat} );
    my $slope_right = ( $p[$right]->{lon} - $p[$iam]->{lon} )
	/ ( $p[$right]->{lat} - $p[$iam]->{lat} );
    my $b_left  = $p[$iam]->{lon} - $slope_left  * $p[$iam]->{lat};
    my $b_right = $p[$iam]->{lon} - $slope_right * $p[$iam]->{lat};

    # finally ready to integrate over left triangle
    my $dlon = ( $lat_ordered[1]->{lon} - $lat_ordered[0]->{lon} ) / 100.0;
    my $dlat = ( $lat_ordered[1]->{lat} - $lat_ordered[0]->{lat} ) / 100.0;

    my $area = 0.;
    for ( my $lat = $lat_ordered[0]->{lat};
	  $lat <= $lat_ordered[1]->{lat}; $lat += $dlat ) {
	my $lon_left  = $slope_left  * $lat + $b_left;
	my $lon_right = $slope_right * $lat + $b_right;
	if ( $lon_left > $lon_right ) {
	    my $t = $lon_left;
	    $lon_left = $lon_right;
	    $lon_right = $t;
	}
	my $coslat = cos( deg2rad( $lat ));
	for ( my $lon = $lon_left; $lon <= $lon_right; $lon += abs($dlon) ) {
	    my $dloncos = $coslat * $dlon;
	    $area += interpolate($shape_p, $lat, $lon)**2
		* deg2rad( abs( $dlat ) )
		* deg2rad( abs( $dloncos ) )
		;
	}
    }

    # now the trapezoid in the middle
    if ( $lat_ordered[1] == $p[$left] ) {
	$left = ($left + 3)%4;
	$slope_left = ( $p[$left]->{lon} - $lat_ordered[1]->{lon} )
	    / ( $p[$left]->{lat} - $lat_ordered[1]->{lat} );
	$b_left = $lat_ordered[1]->{lon}
	- $slope_left * $lat_ordered[1]->{lat};
    }
    else {
	die "confused" unless $lat_ordered[1] == $p[$right];
	$right = ($right + 5)%4;
	$slope_right = ( $p[$right]->{lon} - $lat_ordered[1]->{lon} )
	    / ( $p[$right]->{lat} - $lat_ordered[1]->{lat} );
	$b_right = $lat_ordered[1]->{lon}
	- $slope_right * $lat_ordered[1]->{lat};
    }

    $dlon = ( $lat_ordered[2]->{lon} - $lat_ordered[1]->{lon} ) / 100.0;
    $dlat = ( $lat_ordered[2]->{lat} - $lat_ordered[1]->{lat} ) / 100.0;

    for ( my $lat = $lat_ordered[1]->{lat};
	  $lat <= $lat_ordered[2]->{lat}; $lat += $dlat ) {
	my $lon_left  = $slope_left  * $lat + $b_left;
	my $lon_right = $slope_right * $lat + $b_right;
	if ( $lon_left > $lon_right ) {
	    my $t = $lon_left;
	    $lon_left = $lon_right;
	    $lon_right = $t;
	}
	my $coslat = cos( deg2rad( $lat ));
	for ( my $lon = $lon_left; $lon <= $lon_right; $lon += abs($dlon) ) {
	    my $dloncos = $coslat * $dlon;
	    $area += interpolate($shape_p, $lat, $lon)**2
		* deg2rad( abs( $dlat ) )
		* deg2rad( abs( $dloncos ) )
		;
	}
    }

    # triangle on right.  Integrate lat_ordered[2->3].  That side is
    #  one of the sides
    if ( $lat_ordered[2] == $p[$left] ) {
	$slope_left = ( $lat_ordered[3]->{lon} - $lat_ordered[2]->{lon} )
	    / ( $lat_ordered[3]->{lat} - $lat_ordered[2]->{lat} );
	$b_left = $lat_ordered[3]->{lon}
	- $slope_left * $lat_ordered[3]->{lat};
    }
    else {
	die "confused" unless $lat_ordered[2] == $p[$right];
	$slope_right = ( $lat_ordered[3]->{lon} - $lat_ordered[2]->{lon} )
	    / ( $lat_ordered[3]->{lat} - $lat_ordered[2]->{lat} );
	$b_right = $lat_ordered[3]->{lon}
	- $slope_right * $lat_ordered[3]->{lat};
    }
    $dlon = ( $lat_ordered[3]->{lon} - $lat_ordered[2]->{lon} ) / 100.0;
    $dlat = ( $lat_ordered[3]->{lat} - $lat_ordered[2]->{lat} ) / 100.0;

    for ( my $lat = $lat_ordered[2]->{lat};
	  $lat <= $lat_ordered[3]->{lat}; $lat += $dlat ) {
	my $lon_left  = $slope_left  * $lat + $b_left;
	my $lon_right = $slope_right * $lat + $b_right;
	if ( $lon_left > $lon_right ) {
	    my $t = $lon_left;
	    $lon_left = $lon_right;
	    $lon_right = $t;
	}
	my $coslat = cos( deg2rad( $lat ));
	for ( my $lon = $lon_left; $lon <= $lon_right; $lon += abs($dlon) ) {
	    my $dloncos = $coslat * $dlon;
	    $area += interpolate($shape_p, $lat, $lon)**2
		* deg2rad( abs( $dlat ) )
		* deg2rad( abs( $dloncos ) )
		;
	}
    }
    return $area;
}


sub interpolate
{
    my ($s_p, $lat, $lon) = @_;
    die "bad coordinates ($lat, $lon)\n"
	if $lat < -90
	or $lat > 90
	or $lon < 0
	or $lon > 360;
    my @c = (floor($lat), floor($lat)+1, floor($lon), floor($lon)+1);
    my $slope = retrieve($s_p, $c[1], $c[2]) - retrieve($s_p, $c[0], $c[2]);
    my $a = ($lat - $c[0]) * $slope + retrieve($s_p, $c[0], $c[2]);
    $slope = retrieve($s_p, $c[1], $c[3]) - retrieve($s_p, $c[0], $c[3]);
    my $b = ($lat - $c[0]) * $slope + retrieve($s_p, $c[0], $c[3]);
    $slope = $b - $a;
    return ($lon - $c[2]) * $slope + $a;
}


sub rd_shape_model
    # read evenly distributed list of latitude, longitude, radius
{
    my ($filename) = @_;
    open FILE, '<'.$filename or die "rd_shape_model($filename) open failure\n";
    # use first two lines to determine longitude increment.  These are
    #  redundant w/ 3rd line due to south pole degeneracy
    my ($lat, $lon, $r) = split(' ', <FILE>);
    my $dl = (split(' ', <FILE>))[1] - $lon;
    die "rd_shape_model assumes dl=2 instead of $dl\n" unless $dl == 2;
    my @shape;
    store(\@shape, $lat, $lon, $r);
    store(\@shape, $lat, $lon+$dl, $r);
    while (<FILE>) { store(\@shape, split); }
    close FILE;
    return \@shape;
}


sub retrieve
{
    my ($shape_p, $lat, $lon) = @_;
    return $shape_p->[$lat+90][$lon];
}


sub store
{
    my ($shape_p, $lat, $lon, $r) = @_;
    $shape_p->[int($lat)+90][int($lon)] = $r;
}


my $shape_p = rd_shape_model('areaspud.in');
fill_out_model($shape_p);
print "Give me latitude longitude for each of four corners, corner per line\n";
my @corners;
for (my $c = 0; $c < 4; $c++) {
    my ($lat, $lon) = split(' ', <>);
    my $r = interpolate($shape_p, $lat, $lon);
    $corners[$c] = { lat => $lat, lon => $lon, r => $r };
    print "($lat, $lon) @ $r\n";
}
print integrate($shape_p, \@corners)," km^2 \n";

__END__


# old code fragments
    if (0) {
    
    # now the parallelogram in between the triangles
    $dlon = ( $p[2]->{lon} - $p[1]->{lon} ) / 100.0;
    $dlat = ( $p[2]->{lat} - $p[1]->{lat} ) / 100.0;
    for ( my $lon = $p[1]->{lon}; $lon <= $p[2]->{lon}; $lon += $dlon ) {
	my $lat0 = $slope_1r * $lon + $b_1r;
	my $lat1 = $slope_0q * $lon + $b_0q;
	if ( $lat1 < $lat0 ) {
	    my $t = $lat0;
	    $lat0 = $lat1;
	    $lat1 = $t;
	}
	print "lat0 $lat0 to lat1 $lat1 at lon $lon\n";
	for ( my $lat = $lat0; $lat <= $lat1; $lat += $dlat ) {
	    my $dloncos = cos( deg2rad( $lat ) ) * $dlon;
	    $area += interpolate($shape_p, $lat, $lon)**2 * abs( $dlat )
		* abs( $dloncos )
		* 3.14159 / 180. * 3.14159 / 180.;
#	    print "shape ", interpolate( $shape_p, $lat, $lon ),
#	    " dlat $dlat dloncos $dloncos\n";
	}
    }
    print "area after trapezoid: $area\n";

    # final triangle
    my $a032 = atan2( $p[0]->{lat} - $p[3]->{lat},
		      $p[0]->{lon} - $p[3]->{lon} )
	- atan2(  $p[2]->{lat} - $p[3]->{lat},
		  $p[2]->{lon} - $p[3]->{lon} );
    my $a132 = atan2( $p[1]->{lat} - $p[3]->{lat},
		      $p[1]->{lon} - $p[3]->{lon} )
	- atan2(  $p[2]->{lat} - $p[3]->{lat},
		  $p[2]->{lon} - $p[3]->{lon} );
    print "a032 in radians $a032\n", Math::Trig::rad2rad($a032),"\n";
    print "angles:\n032: ", rad2deg( $a032 ), "\n132: ", rad2deg( $a132 ), "\n";
    exit;
    $dlon = ( $p[3]->{lon} - $p[2]->{lon} ) / 100.0;
    $dlat = ( $p[3]->{lat} - $p[2]->{lat} ) / 100.0;
    for ( my $lon = $p[2]->{lon}; $lon <= $p[3]->{lon}; $lon += $dlon ) {
	my $lat0 = $slope_23 * $lon + $b_23;
	my $lat1 = $slope_0q * $lon + $b_0q;
	if ( $lat1 < $lat0 ) {
	    my $t = $lat0;
	    $lat0 = $lat1;
	    $lat1 = $t;
	}
	print "lat0 $lat0 to lat1 $lat1 at lon $lon\n";
	for ( my $lat = $lat0; $lat <= $lat1; $lat += $dlat ) {
	    my $dloncos = cos( deg2rad( $lat ) ) * $dlon;
	    $area += interpolate($shape_p, $lat, $lon)**2 * abs( $dlat )
		* abs( $dloncos )
		* 3.14159 / 180. * 3.14159 / 180.;
#	    print "shape ", interpolate( $shape_p, $lat, $lon ),
#	    " dlat $dlat dloncos $dloncos\n";
	}
    }
}


sub new
{
    my ( $shape_p, $corners_p ) = @_;
    my @p = sort { $a->{lon} <=> $b->{lon} } @$corners_p;

    # left-most triangle.
    #  first, is angle 102 or 103 larger?
    my $a102 = atan2( $p[2]->{lat} - $p[1]->{lat},
		      $p[2]->{lat} - $p[1]->{lat} );
    my $a103 = atan2( $p[3]->{lat} - $p[1]->{lat},
		      $p[3]->{lat} - $p[1]->{lat} );
    print "angles:\n102: ", rad2deg( $a102 ), "\n103: ", rad2deg( $a103 ), "\n";
    my ( $q, $r );
    if ( abs( $a102 ) > abs( $a103 ) ) { $q = $p[2]; $r = $p[3]; }
    else { $q = $p[3]; $r = $p[2]; }

    #  parameterize those sides y=mx+b:
    my $slope_01 = ( $p[1]->{lat} - $p[0]->{lat} )
	/ ( $p[1]->{lon} - $p[0]->{lon} );
    my $b_01 = $p[0]->{lat} - $slope_01 * $p[0]->{lon};
    my $slope_0q = ( $q->{lat} - $p[0]->{lat} )
	/ ( $q->{lon} - $p[0]->{lon} );
    my $b_0q = $p[0]->{lat} - $slope_0q * $p[0]->{lon};
    my $slope_1r = ( $r->{lat} - $p[1]->{lat} )
	/ ( $r->{lon} - $p[1]->{lon} );
    my $b_1r = $p[1]->{lat} - $slope_1r * $p[1]->{lon};

    my $dlon = ( $p[1]->{lon} - $p[0]->{lon} ) / 100.0;
    my $dlat = ( $p[1]->{lat} - $p[0]->{lat} ) / 100.0;

    my $area = 0.;
    for ( my $lon = $p[0]->{lon}; $lon <= $p[1]->{lon}; $lon += $dlon ) {
	my $lat0 = $slope_01 * $lon + $b_01;
	my $lat1 = $slope_0q * $lon + $b_0q;
	if ( $lat1 < $lat0 ) {
	    my $t = $lat0;
	    $lat0 = $lat1;
	    $lat1 = $t;
	}
	print "lat0 $lat0 to lat1 $lat1 at lon $lon\n";
	for ( my $lat = $lat0; $lat <= $lat1; $lat += $dlat ) {
	    my $dloncos = cos( deg2rad( $lat ) ) * $dlon;
	    $area += interpolate($shape_p, $lat, $lon)**2 * abs( $dlat )
		* abs( $dloncos )
		* 3.14159 / 180. * 3.14159 / 180.;
#	    print "shape ", interpolate( $shape_p, $lat, $lon ),
#	    " dlat $dlat dloncos $dloncos\n";
	}
    }
    print "triangle1: $area\n";


    print "final area: $area\n";

    warn "not tested yet!\n";
    return $area;
}
